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5  INTRODUCTION 

5.1  Nature  of  the  problem 

The  potential  benefits  of  early  diagnosis  of  cancer  were  recognized  many  years  ago,  before 
soft  tissue  imaging  was  available.  This  goal  was  behind  the  first  efforts  to  apply  ultrasound 
to  the  problem  of  the  detection/diagnosis  of  breast  cancer.  Since  then,  many  investigators 
have  devoted  effort  to  this  problem  because  of  the  known  advantages  of  ultrasound:  it  is 
non-ionizing,  relatively  inexpensive  and  uses  widely  available,  portable,  equipment.  Today, 
with  the  success  of  x-ray  mammography  as  an  early  screening  tool,  there  is  still  room  for 
improved  methods,  since  there  is  disagreement  regarding  the  use  of  ionizing  radiation  for 
screening,  or  routine  exams.  Also,  the  fibrous  tissue  of  dense  breasts  gives  poor  results  in 
conventional  mammograms. 

Ultrasound  is  now  available  as  an  adjunctive  modality  in  many  breast  clinics,  where 
it  is  used  to  determine  if  masses  with  smooth  borders  are  cystic  or  solid  and  to  examine 
dense  young  breasts.  Biopsy  is  still  used  often  to  determine  malignancy.  The  cystic/solid 
determination  can  be  a  problem  in  borderline  cases,  because  of  the  tendency  of  some 
ultrasonic  systems  to  “fill  in”  the  echo  free  space  with  artifactual  echoes  and  thus  make 
a  cyst  appear  to  be  solid.  There  is  a  need  for  better  ultrasound  systems,  as  well  as  to 
have  a  more  general-purpose  imaging  modality  available  than  the  x-ray,  particularly  if  this 
modality  had  the  advantages  listed  for  ultrasound. 

Several  factors  have  been  identified  in  the  literature  as  those  contributing  to  poor  res¬ 
olution  of  ultrasound  images.  In  B-Scan  images,  finite  signal  bandwidth  of  the  ultrasonic 
transducers  is  a  major  reason  for  low  resolution  in  the  temporal  axis,  whereas  the  non 
negligible  beam  width  highly  contributes  to  that  in  the  lateral  direction  [18],  [7],  [14],  [44], 
[21].  The  resolution  also  depends  on  the  frequency  at  which  the  imaging  system  operates. 
In  addition  to  equipment  limitations,  there  are  factors  originating  from  the  nature  of  the 
tissue  being  imaged.  Phase  aberrations  and  velocity  variations  arising  from  acoustic  in¬ 
homogeneity  of  tissues  are  two  of  the  important  causes,  not  only  for  low  resolution  but 
also  for  low  contrast  in  images  [39],  [11],  [6].  The  observed  ultrasonic  image  can,  therefore, 
be  considered  as  a  distorted  version  of  the  true  tissue  image,  where  the  axial  distortion  is 
dominated  by  the  pulse-echo  wavelet  of  the  imaging  system  and  the  lateral  distortion  by 
the  lateral  beam  profile. 
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5.2  Background  of  previous  work 

The  problem  of  compensating  for  imaging  distortions  has  been  very  active  research  area, 
underscoring  the  importance  in  improving  diagnostic  quality  of  ultrasonic  images.  There 
have  been  several  approaches  to  cancel  out  the  aberration  effects,  but  there  is  no  consensus 
as  to  the  best  way  to  achieve  it  [10],  [26],  [25],  [38].  In  [10]  it  was  proposed  estimating 
the  differences  in  arrival  times  between  two  adjacent  receiving  element  locations  using 
cross  correlation  techniques;  these  results  were  used  to  modify  phasing  characteristics  of 
the  transducers  for  ensuing  scans.  In  [25],  using  an  idea  adapted  from  optics,  phasing 
characteristics  were  determined  using  speckle  brightness  as  a  measure  of  image  quality. 

The  availability  of  convolution  models,  such  as  in  [17],  [9],  for  ultrasound  image  for¬ 
mation  and  the  wide  availability  of  digital  computers  has  given  an  added  importance  to 
discrete-time  deconvolution  methods,  as  a  means  of  improving  images  beyond  the  capabil¬ 
ities  of  hardware.  A  few  researchers  have  investigated  the  true  2-D  deconvolution  of  RF 
images  [9],  [18],  [32],  whereas  most  of  the  published  works  are  on  1-D  techniques  [44],  [21], 
123],  [19],  [14], 

Since  the  resolution  along  the  lateral  direction  was  much  worse  than  that  along  the 
temporal  axis,  a  number  of  attempts  were  focussed  on  deconvolution  of  lateral  image  lines. 
In  [44],  a  B-mode  image  was  considered  to  be  an  ensemble  of  lateral  lines  corresponding 
to  lateral  slices  through  the  envelope  detected  image,  at  given  times  (depth).  Observing 
that  the  point-spread  function  of  a  typical  pulse-echo  imaging  system  is  highly  elongated 
along  the  lateral  direction,  they  hypothesized  that  lateral  image  lines  can  be  approximately 
described  by  a  1-D  convolution  model.  Their  model  consisted  of  two  1-D  terms:  a  signal  of 
interest  called  the  tissue  reflectance  and  a  blurring  kernel  in  the  lateral  direction  called  the 
lateral  point  spread  function.  The  latter  function  was  defined  to  be  the  laterally  varying 
component  of  the  2-D  point  spread  function,  whose  axial  variation  had  been  approximated 
by  a  Dirac  delta  function.  The  problem  of  resolution  enhancement  was  posed  as  one  of 
extracting  the  the  tissue  reflectance  from  the  observed  image,  assuming  a  perfect  knowledge 
of  the  lateral  point  spread  function.  Using  a  Gaussian  shaped  hypothetical  lateral  point 
spread  function,  it  was  shown  that  at  the  best  signal  to  noise  ratio  that  can  be  expected 
from  ultrasound  images,  deconvolution  will  lead  to  a  resolution  enhancement  of  no  better 
than  2.0.  The  definition  of  the  resolution  was  based  on  the  reciprocal  of  the  effective  width 
of  the  lateral  point  spread  function.  However,  the  amount  of  improvement  was  also  reported 
to  be  dependent  on  the  exact  shape  of  the  lateral  point  spread  function.  These  figures  were 
found  to  be  in  agreement  with  the  empirical  numbers  reported  in  [14].  Several  others  have 
reported  results  on  1-D  lateral  deconvolution  [43],  [33],  [21]  ,  where  in  [33]  it  was  concluded 
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that  the  computational  effort  on  lateral  deconvolution  was  wasted  because  of  the  very  low 
resolution  enhancement  they  could  obtain  at  the  expense  of  introducing  more  artifacts. 

The  drawbacks  of  the  lateral  deconvolution  techniques  discussed  above  are  the  following. 
Although  the  quantity  displayed  on  an  ultrasound  imager  is  the  envelope  of  the  received  RF 
signal,  the  image  formation  process  actually  occurs  in  the  RF  domain.  In  [32]  it  is  pointed 
out  that  in  general  the  convolutional  model  of  image  formation  in  the  RF  domain  does 
not  hold  for  envelope  detected  signals.  It  is  concluded  that  the  1-D  lateral  deconvolution 
on  envelope  succeeds  only  in  the  special  case  where  no  phase  interference  from  nearby 
reflectors  is  present. 

All  the  papers,  except  [21],  on  lateral  deconvolution  mentioned  above,  relied  on  mea¬ 
suring  the  lateral  beam  profile  to  be  used  in  computations.  To  obtain  the  true  tissue  image 
from  the  distorted  observation,  one  requires  quantitative  information  on  the  complex  beam 
shape  and  acoustic  velocity  variations  in-vivo,  which  are  impossible  to  measure  directly.  In 
highly  simplified  situations  such  as  measurements  of  wire  targets  under  water,  it  is  possible 
to  get  that  information  reliably.  However,  in  the  case  of  in-vivo  tissue  targets,  such  infor¬ 
mation  is  generally  unavailable.  The  measurements  done  under  water  are  not  valid  with 
clinical  images,  even  when  the  imaging  system  used  is  the  same,  because  of  the  effects  of 
phase  aberrations,  nonlinearities  and  dispersive  attenuation  introduced  by  the  tissues  [21], 
[11],  [39]  ,  [6].  thus  dramatically  limiting  the  clinical  applicability  of  those  methods. 

In  [21],  a  line-by-line  lateral  deconvolution  technique  which  does  not  require  the  complex 
beam  shape  in  tissue  or  phase  information  on  adjacent  lateral  lines,  was  proposed.  This 
method  too,  however,  worked  on  amplitude  detected  B-scan  images.  It  hinges  on  the  key 
assumption  that  the  transfer  function  of  the  imaging  system  along  a  lateral  image  line 
can  be  approximated  by  the  smoothed  Fourier  transform  of  the  lateral  image  line  itself. 
A  convolution  relationship  between  the  envelope  of  a  lateral  point  spread  function  and  a 
slowly  varying  envelope  of  the  tissue  response  has  been  tacitly  assumed.  Thus  this  method 
is  subject  to  all  the  limitations  implied  by  the  above  assumptions. 

In  spite  of  the  fact  that  temporal  resolution  in  an  ultrasound  image  is  much  higher 
than  the  lateral  resolution,  axial  deconvolution  is  still  of  importance.  Besides  the  obvious 
advantage  in  improved  axial  resolution,  the  removal  of  the  effect  of  the  ultrasound  pulse 
echo  wavelet  (through  axial  deconvolution)  will  tend  to  make  the  appearance  of  images  more 
uniform  over  different  subjects  [16],  thus  simplifying  the  diagnosis  procedures.  As  the  shape 
of  the  pulse  echo  wavelet  changes  with  propagation  due  to  dispersive  attenuation,  a  first 
step  in  axial  deconvolution  often  involves  the  estimation  of  the  pulse  in  tissue.  Parametric 
modeling  of  speckle-only  image  lines  has  been  proposed  [37],  [20],  [16].  However,  such 
approaches  are  limited  by  the  problems  such  as  the  model  order  selection,  associated  with 
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parametric  modeling.  In  [19]  a  Kalman  filter  technique  was  applied  to  estimate  pulse  echo 
wavelets  as  well  as  to  simultaneously  improve  the  axial  resolution.  The  success  of  the 
method  was  reported  to  be  dependent  on  the  SNR  of  the  observations  and  the  accuracy  at 
which  the  observations  could  be  modeled.  A  non-parametric  approach  for  the  estimation 
of  the  pulse  was  proposed  in  [15],  where  the  minimum-phase  equivalent  of  the  pulse-echo 
wavelet  was  separated  from  the  tissue  response.  However,  quite  often  pulse-echo  wavelets 
and  lateral  kernels  are  non- minimum  phase  signals,  thus  limiting  the  generality  of  this 
approach. 

5.3  Methods  of  approach 

In  the  past  yeax  we  introduced  a  novel  non-parametric  framework  for  deconvolution  of 
B-scan  images  [1],  [2],  [3],  [4],  [5].  We  first  developed  a  model  for  the  rf  image,  and  then 
reconstructed  distortions  using  higher-order  statistics  of  the  measured  image  lines.  Based 
on  the  estimated  distortions  we  performed  deconvolution  of  the  corresponding  images  and 
demonstrated  that  the  resolution  of  ultrasound  images  of  tissue  mimicking  phantom  as  well 
as  human  tissue  images  was  significantly  improved.  In  the  past,  estimation  of  distortions 
has  been  carried  out  using  exclusively  second  order  statistics  (autocorrelation).  Autocor¬ 
relation,  however,  can  recover  only  the  minimum  phase  quivalent  of  the  true  distortions, 
because  it  is  blind  to  phase.  We  also  estimated  distortions  using  second  order  statistics 
of  the  image  lines.  We  showed  that  although  these  estimates  did  lead  to  resolution  im¬ 
provement,  the  amount  of  improvement  was  less  significant  that  the  one  obtained  with  the 
higher-order  statistics  based  estimates. 

Image  formation  process  in  the  RF-domain  is  described  by  a  2-D  convolutional  model, 
where  the  attenuation  of  the  pulse-echo  wavelet  and  beam  aberration  effects  can  be  indi¬ 
rectly  incorporated  [17],  [18],  [16].  Two  1-D  blurring  kernels,  corresponding  to  axial  and 
lateral  directions,  is  hypothesized  to  represent  distortions  along  respective  axes.  The  axial 
distortion  kernel  includes  the  blurring  effects  due  to  the  finite  bandwidth  of  the  transducer, 
and  dispersive  attenuation  of  the  pulse-echo  wavelet  in  tissue.  The  lateral  distortion  kernel 
represents  the  convolutional  components  of  lateral  blurring  due  to  the  complex  beam  pat¬ 
terns.  Formalizing  a  definition  for  resolution,  we  show  that  compensation  for  the  effects  of 
the  blurring  kernels  improves  the  resolution  of  the  image. 

The  proposed  method  has  the  advantage  of  being  able  to  estimate  these  kernels  at 
each  image  line,  axial  and  lateral,  thus  capturing  the  variations  within  the  image.  Since 
the  estimations  are  based  on  higher-order  statistics  [24]  of  RF-data,  the  estimated  kernels 
are  robust  to  additive  observation  noise  and  also  have  correct  phase.  To  the  best  of  our 
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knowledge,  the  method  we  proposed  [1]  is  the  first  one  to  in-vivo  estimate  the  distortion 
kernels  with  their  true  phase,  as  opposed  to  conventional  methods  that  estimate  minimum- 
phase  equivalent  of  kernels. 

6  BODY 

6.1  Modeling  the  rf  image 

During  an  ultrasonic  investigation,  a  three-dimensional  pulsed  pressure  field  is  emitted 
into  the  tissue.  The  field  interacts  with  the  tissue  and  part  of  it  is  reflected,  scattered,  and 
subsequently  received  by  the  transducer.  Under  the  assumptions  of  linear  propagation  and 
weak  scattering,  an  expression  for  the  received  pressure  field  was  derived  in  [17],  using  the 
first  order  Born  approximation.  Absorption  effects  were  neglected.  The  equation  has  been 
expressed  as  a  convolutional  model  in  the  following  form: 

y{r2,t)  =  Vp^{t)  *t  /(ri)  *r  hpe{ri,r2,t)  -\-  w{r2,t),  (1) 

where: 

•  ri,  r2  are  vectors  denoting  the  location  of  the  scatterer  and  the  transducer,  respectively; 

•  “  *t  ”  and  denote  time  and  spatial  convolution,  respectively; 

•  /(ri)  originates  from  the  inhomogeneities  in  the  tissue  due  to  density  and  propagation 
velocity  perturbations  above  their  mean  levels,  giving  rise  to  the  back  scattered  signal  {tis¬ 
sue  response).] 

•  Vpe{t)  is  the  pulse-echo  wavelet  that  accounts  for  the  transducer  excitation  and  the  im¬ 
pulse  responses  during  emission  and  reception  of  the  pulse; 

•  hpe{ri,r2,  t)is  the  modified  pulse-echo  spatial  impulse  response  that  relates  the  transducer 
geometry  and  the  spatial  extend  of  the  scattered  field.  The  computation  of  hpe(ri,r2,t) 
is  based  on  the  approach  described  in  [41],  [35].  Convolutional  components  of  aberrations 
and  dispersive  attenuation,  which  introduce  spatially  varying  effects  to  the  process,  may 
be  incorporated  in  this  already  spatially  varying  kernel. 

•  w{r2,t)  represents  measurement  noise  and  the  unmodeled  dynamics  of  the  image  forma¬ 
tion  process. 

The  problem  of  extracting  the  tissue  response  /(ri)  from  the  observation  y{r2,t)  is  a 
deconvolution  problem.  Since  a  B-mode  image  is  a  mapping  from  the  3-D  tissue  space  to 
the  2-D  space  of  the  display,  the  solution  is  not  unique  in  general.  This  non-uniqueness 
should  be  obviated  by  making  reasonable  assumptions  about  the  3-D  structure  of  the  tissues 
being  imaged  [9].  An  assumption  implicit  in  cases  where  the  deconvolution  is  done  using 
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kernels  confined  to  the  imaging  plane  is  that  all  image  features  in  the  imaging-plane  extend 
perpendicular  to  the  corresponding  plane  in  the  tissue  space,  so  as  to  make  the  height  of 
extension  the  effective  height  of  the  beam. 

The  convolution  model  of  (19)  expresses  the  fact  that  the  received  signal  at  the  trans¬ 
ducer  site  r2  is  a  result  of  linear  spatio-temporal  interaction  between  the  signal  of  interest 
/(ri)  and  a  distortion  kernel.  Thus,  the  measured  signal  contains  a  distorted  version  of 
the  true  tissue  response  /(ri ) .  Deconvolving  these  kernels  should  improve  image  resolution 
and  contrast.  It  should  also  remove  aberration-induced  artifacts  that  result  from  changing 
beam  profiles  inside  the  tissue. 

As  discussed  in  the  introduction,  efforts  to  carry  out  this  deconvolution  have  been 
hampered  by  the  difficulties  in  measuring  the  modified  spatial  impulse  response  of  the 
imaging  system.  Underwater  measurements  using  simplified  targets  would  not  reveal  any 
significant  changes  undergone  by  the  interrogating  beam  in  tissue. 

Our  goal  here  is  actually  to  identify  the  combination  of  Vpe{t)  and  the  spatially  varying 
kernel  /ipe(ri,r2,t),  and  subsequently  cancel  it  from  the  image  in  order  to  improve  lateral 
as  well  as  axial  resolution.  Let  us  combine  both  smoothing  kernels  Vpe{t)  and  hpe(ri,r2it) 
in  one  spatially  and  temporally  varying  kernel,  /j(ri,r2,t),  which  we  are  going  to  refer  to 
as  the  ultrasonic  system  impulse  response.  For  discrete  time,  and  for  some  fixed  transducer 
location  (19)  is  equivalent  to  the  following  two-dimensional  convolutional  model 

yihn)  =  +  (2) 

i  3 

where  y(/,  n)  represents  the  sample  from  the  /— th  A-line  at  discrete  time  n.  The  goal  here 
is  to  identify  the  time  varying  ultrasonic  system  response  h(l,n),  and  recover  the  tissue 
response,  f{i,j),  from  the  noisy  measurement  y{k,l). 

However,  we  will  not  attempt  a  true  2-D  deconvolution  in  this  paper.  As  is  commonly 
done  in  ultrasound  deconvolution  literature  [15],  [23],  [18],  [19],  we  assume  that  an  RF 
A-line  can  be  expressed  as  a  convolution  between  two  1-D  axial  terms:  a  hypothetical 
tissue  response  and  a  distortion  kernel.  This  view  is  not  unique  to  the  image  deconvolution 
literature;  ultrasound  Doppler  systems  and  tissue  attenuation  estimation  techniques  tacitly 
depend  on  it  [15].  Reducing  the  problem  to  a  1-D  deconvolution  is  analogous  to  the  original 
decomposition  of  the  true  3-D  deconvolution  problem  in  to  a  2-D  one. 
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6.1.1  Assuming  white  tissue  response:  In  vivo  distortion  estimation  using 
higher-order  statistics  (  HOS) 

In  the  following  we  will  treat  each  image  line  (either  in  ateral  or  axial  direction)  separately 
assuming  the  1-D  model 


yi{n)  =  hi{n)  *  fi{n)  +  Wi{n),  i  =  1,2, ...,  (3) 

where  i  is  the  A-line  index;  fi{n)  is  the  axial  or  lateral  tissue  response;  hi{n)  is  the  axial  or 
lateral  distortion  kernel;  which  describes  the  distortion  associated  with  the  i-th  line;  and  n 
denotes  discrete  time.  We  assumed  that: 

(Al)  hi(n)  is  deterministic,  possibly  non-minimum  phase, 

(A2)  fi{n)  is  stationary,  white,  independent  identically  distributed  (i.i.d.),  zero-mean  non- 
Gaussian, 

(A3)  Wi{n)  is  white  zero-mean  Gaussian,  and  independent  of  fi{n). 

If  we  consider  a  region  of  the  image  that  contains  speckle  only,  these  assumptions  should 
hold  reasonably  well.  In  an  image  that  contains  inhomogeneities,  however,  the  whiteness 
assumption  about  /t(n)  will  not  be  valid. 

The  third  order  cumulant  of  the  observation  j/,(n)  is  defined  as  [24]  : 

%  P)  =  E{yi{n)yi{n  -f  T)i/;(n  +  p)}  (4) 

The  bispectrum  of  the  yi{n)  is  defined  to  be  the  Fourier  transform  of  the  third  order 
cumulant.  Under  assumptions  (Al)  —  (A3),  the  bispectrum  of  y,(n)  is  given  by  [24]: 


Cy-{L0i,uj2)  —  C -\-0j2)  -b  C'u,;(a;i,a;2)  (5) 

where  (7/;(a;i,u>2)  is  the  bispectrum  of  /i(n),  Hi{u))  is  the  spectrum  of  hi{n)  and  C.w^{u!i,otJ2) 
is  the  noise  bispectrum.  If  the  additive  noise  is  zero-mean  Gaussian,  then  Cuji(toi,u}2)  —  0 
[24]  and  (4)  becomes, 

Cy-(uJl,LL)2)  =  C fX^l,L02)Hi{L0i)Hi{u)2)H*(uii  +  ^2)  (6) 

The  bicepstrum  is  defined  as  the  cepstrum  of  the  bispectrum,  from  which  we  get: 

6y<(mi,m2)  =  6/i(mi,m2)  +  Cft.(mi,m2),  (7) 


where  bf-{mi,m2),  c/j;(mi,m2)  are  the  bicepstra  of  /j(n)  and  hi{n),  respectively. 
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If  /j(n)  is  stationary  independent  identically  distributed  (i.i.d.),  its  third  order  spectrum 
is  flat  and  equal  to  the  skewness,  7/;,  of  the  process.  Therefore,  its  bicepstrum  will  be  an 
impulse  located  at  the  origin.  In  that  case,  using  bicepstral  values  along  the  main  axes 
except  at  the  origin,  we  can  reconstruct  a  scaled  and  shifted  version  of  hi(n)  as: 


hi{n)  = 


(8) 


where 


Cfti(m) 


by-{m,0)  m  >  0 

0  m  =  0 


(9) 


I  rn<d. 

We  assume  that  a  similar  1-D  model  holds  in  the  lateral  direction  of  the  RF-image. 
Then,  a  similar  procedure  can  be  followed  to  estimate  lateral  distortion  kernels  at  each  line. 
The  kernels  thus  estimated  will  include  the  convolutional  components  of  aberration  as  well. 
Therefore,  axial  deconvolution  followed  by  lateral  deconvolution,  or  various  combinations 
thereof,  should  give  us  distortion  compensated  RF  images  which  will  have  higher  resolution 
and  contrast.  The  ability  to  estimate  and  remove  beam  distortion  effects  (due  to  aberration) 
is  seen  as  a  major  advantage  of  this  method  over  other  non-parametric  techniques  [15],  [21]. 
The  method  proposed  in  [15],  as  is,  can  not  estimate  non- minimum  phase  signals,  while 
that  of  [21]  has  been  designed  for  envelope  detected  signals. 


6.1.2  Assuming  non-white  tissue  response:  In  vivo  distortion  estimation  using 
HOS 

As  it  was  already  mentioned  the  tissue  response  fi{n)  will  be  non- white  if  variations  in 
tissue  structure  occur.  A  detailed  study  of  this  case  can  be  found  in  [2]  (see  also  Appendix 
A).  The  most  common  way  to  model  a  non- white  process  is  consider  it  as  the  convolution 
of  a  white  noise  term  with  a  deterministic  kernel.  This  kernel  is  referred  to  as  the  “color” 
of  the  process.  It  describes  the  correlation  between  the  samples  of  the  process,  and  can  be 
used  as  an  indication  of  the  the  deviation  of  the  process  from  whiteness. 

If  the  process  described  in  the  previous  Section  is  applied  when  the  tissue  response 
is  nonwhite,  the  distortions  reconstructed  from  each  image  line  will  contain  as  convolu¬ 
tional  component  the  color  of  the  tissue  response.  These  distortions  will  be  referred  to  as 
“combined  kernels”. 

Combined  distortion  kernels  ki{n)  and  kj{n)  obtained  at  two  different  but  closely  spaced 
lines  “i”  and  can  be  written  as: 


ki[n)  =  h{n)  *  ti(n) 


(10) 
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kj{n)  =  h{n)  *  tj{n)  (11) 

where  h{n)  denotes  the  true  distortion  kernel,  which  is  assumed  to  be  invariant  between 
the  two  locations  “i”  and  “j”;  ti{n)  and  tj{n)  respectively  represent  the  color  of  the  tissue 
responses  at  “i”  and  “j”. 

The  true  distortion  h{n)  and  the  colors  tj{n)  can  be  computed  via  the  blind 

deconvolution  method  presented  in  [31]. 


6.1.3  Assuming  white  tissue  response:  In  vivo  estimation  of  the  minimum 
phase  equivalent  of  distortion  using  second-order  statistics  (SOS) 

Using  the  model  of  eq.  (3),  second  order  statistics  were  also  used  to  estimate  imaging 
distortions  [12],  [13]. 

Transforming  (3)  in  the  autocorrelation  domain  (second-order  statistic)  we  get: 

ryti'r)  =  E{yi{n)yi(n  +  r)}  =  7/'  ^  hi{n)hi{n  -f  r)  +  7r<^(r),  (12) 

n 

where  is  the  variance  of  /j(n),  7^'  is  the  variance  of  the  noise,  and  ^(r)  is  the  unit 
impulse.  The  power  cepstrum  is  defined  as  the  inverse  Fourier  transform  of  the  logarithm 
of  the  Fourier  transform  of  the  autocorrelation.  Assuming  that  the  noise  level  is  low  enough, 
in  the  power  cepstrum  domain  we  get: 

yi{k)  ~  hi{k)^  ^  7^  0,  (13) 


where  hi{k)  is  the  power  cepstrum  of  the  distortion  kernel  associated  with  the  th  line. 
From  its  power  cepstrum,  we  can  reconstruct  the  minimum  phase  equivalent  of  hi{n)  [27] 
as: 

hi{n)  =  F~'^[exp{F{yi{k)w{k)]}],  (14) 

where  F].]  and  F’“^[.]  denote  forward  and  inverse  Fourier  transforms,  and 


w{k)  = 


1,  k>0 
0,  otherwise. 


(15) 


6.2  Experiments 

The  goals  of  our  experimental  work  were  to: 

1.  Compare  the  axial  distortions,  estimated  based  on  higher-order  statistics  (HOS-based 
distort  ions),  with  experimentally  obtained  ones.  Validate  the  imposeed  assumptions. 
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2.  Estimate  HOS-based  distortions  from  clinical  ultrasound  images.  Deconvolve  ultra¬ 
sound  images  using  the  HOS-based  estimated  distortions  based  on  the  whiteness 
assumptions  about  the  tissue  response  (see  Section  6.1.1).  Quantify  resolution  im¬ 
provement. 

3.  Estimate  HOS-based  distortions  from  clinical  ultrasound  images.  Deconvolve  ultra¬ 
sound  images  using  the  HOS-based  estimated  distortions  based  on  the  non-whiteness 
assumptions  about  the  tissue  response  (see  Section  6.1.2).  Quantify  resolution  im¬ 
provement.  Determine  if  there  is  any  difference  in  the  final  deconvolution  result 
depending  on  the  assumption  about  the  whiteness  and  non-whiteness  of  the  tissue 
response. 

4.  Deconvolve  clinical  images  using  distortions  estimated  based  on  second-order  statis¬ 
tics  (SOS-based  distortions)  as  outlined  in  Section  6.1.3).  Quantify  resolution  im¬ 
provement. 

5.  Compare  resolution  improvement  of  deconvolution  with  HOS-based  and  SOS-based 
distortions. 

Towards  the  above  goals  we  collected  the  following  data: 

•  Water-tank  measurements: 

Experiment  A 

We  obtained  B-Scan  images  of  an  ATS  model  532  contrast  resolution  phantom 
which  was  positioned  in  a  water  tank.  The  target  area  consisted  of  the  tissue 
mimicking  background  which  had  a  scatterer  density  of  32  scatterers/mm^.  The 
transducer  we  used  was  a  model  GEl  04682  curved  device  with  a  nominal  center 
frequency  of  3.5  MHz.  The  nominal  focal  zone  of  the  transducer  was  6-13  cm.  A 
stepper  motor  with  a  step  size  0.025mm  controlled  the  position  of  the  transducer 
under  the  guidance  of  a  personal  computer.  Data  acquisition  was  done  through 
a  LeCroy  model  9450A  dual  channel  digital  oscilloscope  connected  to  the  PC  by 
a  GPIB  interface.  Transducer  was  moved  across  the  scanning  plane  in  steps  of 
0.25mm,  and  RE  echos  were  sampled  at  a  13.3  MHz  rate. 

Experiment  B 

A  long  piece  of  freshly  peeled  wire  of  diameter  0.812mm  was  placed  underwater 
so  that  it  was  parallel  to  the  transducer  surface  and  perpendicular  to  the  scan¬ 
ning  plane.  It  was  kept  inside  the  focal  zone  of  the  transducer,  at  a  distance 
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8cm  away  from  the  transducer  surface.  The  experimental  setup  used  here  is 
the  same  as  in  experiment  (A).  The  wavelet  reflected  from  the  wire  surface  was 
recorded  and  taken  to  be  an  approximate  estimate  of  the  pulse-echo  wavelet  of 
the  imaging  system. 

•  Measurements  from  Clinical  Equipment: 

Experiment  C 

To  demonstrate  the  performance  of  our  method  on  data  from  more  realistic 
equipment,  we  imaged  the  same  ATS532  phantom  using  a  linear  array  sector 
scan  transducer  on  a  model  UltraMark-9  clinical  imaging  system  manufactured 
by  Advanced  Technology  Laboratories,  Seattle,  U.S.A.  The  target  area  consisted 
of  two  cylinders  with  scatterer  densities  4  and  8  scatterers/mm^,  embedded 
in  a  tissue  mimicking  background  of  32  scatterers/mm^.  The  nominal  center 
frequency  of  the  scanner  was  3.5  MHz;  the  field  of  vision  was  60°.  Data  were 
sampled  at  a  rate  of  12  MHz.  No  TGC  was  employed. 

Experiment  D 

To  demonstrate  the  performance  of  our  method  under  clinical  conditions,  we 
obtained  liver  scans  of  patients  imaged  by  ultrasonologists  at  the  Thomas  Jef¬ 
ferson  University  Hospital,  Philadelphia,  on  the  same  imaging  system  described 
under  experiment  (C).  The  patients  have  been  diagnosed  with  hypoechoic  mul¬ 
tiple  liver  metastatic  tumors.  No  T.G.C.  had  been  applied;  focus  of  the  imaging 
system  at  transmission  was  set  at  2cm,  and  dynamic  focusing  was  employed 
with  the  receiving  mode. 

Experiment  E 

Clinical  breast  data  consists  of  ultrasound  images  acquired  by  radiologists,  from 
patients  referred  to  the  Department  of  Radiology,  Thomas  Jefferson  Universi  ty 
Hospital,  Philadelphia  19107.  The  imaging  system  was  an  Ultramark  9  HDl 
system  with  RF  data  acquisition  capabilities,  provided  to  Thomas  Jefferson 
University  Hospital  by  Advanced  Technology  Laboratories.  A  flat,  linear  array 
transducer  of  center  frequency  7.5  MHz  was  used  to  collect  data  (RF-echoes). 
The  sampling  rate  of  the  data  acquisition  system  was  at  20  MHz.  Five  breast 
images  (denoted  as  Patient  1-  Patient  5)  with  tumors  in  them  are  considered 
here. 
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6.3  Results 

6.3.1  Water-tank  experiments  and  the  accuracy  of  our  estimation  method 

In  order  to  demonstrate  the  validity  of  our  kernel  estimation  procedure,  we  used  the  data 
from  underwater  measurements,  i.e.  experiments  (A)  and  (B).  A  part  of  the  RF  image  data 
that  we  gathered  using  the  circular  transducer  in  experiment  (A)  is  shown  in  Fig.  1(a), 
where  the  logarithm  of  the  envelope  has  been  used  for  display  purposes. 

Under  the  assumption  that  the  axial  blurring  kernel  of  an  image  line  does  not  signifi¬ 
cantly  depend  on  the  lateral  location  of  the  line,  data  from  several  nearby  axial  RF  lines 
can  be  used  to  make  a  longer  data  vector,  which  will  enable  us  to  obtain  better  cumulant 
estimates.  In  order  to  minimize  the  effect  of  attenuation  on  our  estimations,  from  each 
A-line  i  we  considered  data  segments  yi{k)  of  length  not  more  than  2N  samples.  The 
number  of  adjacent  line  segments  yi{k)  concatenated  was  in  the  range  i  =  1  —  10.  In  the 
axial  kernel  estimations,  we  used  M  =  10  and  N  =  64. 

Axial  blurring  kernels  estimated  with  the  method  of  Section  6.1.1  from  different  regions 
of  the  image  in  Fig  1(a)  are  shown  in  Fig.  2(a),  in  dotted  lines.  All  of  these  kernels  have 
been  estimated  from  axial  data  obtained  at  the  same  depths  of  the  image,  but  at  different 
lateral  locations.  The  mean  axial  blurring  kernel,  Vm{t)^  which  was  computed  as  the  average 
of  estimated  kernels,  is  indicated  by  the  solid  line  while  the  measured  pulse  echo  wavelet 
is  indicated  by  the  dashed  line. 

It  can  be  seen  that  all  of  the  estimated  kernels  possess  a  similar  structure,  resembling 
a  typical  ultrasound  pulse-echo  wavelet.  This  is  to  be  expected  since  the  component  due 
to  the  pulse-echo- wavelet  (upe(t))  dominates  the  axial  blurring  kernel  [18].  The  variation 
among  the  estimated  kernels  may  be  attributed  to  the  statistical  estimation  errors,  de¬ 
viation  of  the  scatterer  response  from  a  statistically  white  response  (which  violates  our 
assumption  (A2))  and,  contributions  from  the  effects  of  the  medium  such  as  aberration 
and  dispersive  attenuation  which  have  spatially  varying  characteristics. 

There  is  a  reasonable  match  between  the  measured  and  average  estimated  kernels  (Fig. 
2),  within  the  limits  of  the  accuracy  of  the  experimentation.  Looking  at  their  spectra  (see 
Fig.  2(b))  we  can  clearly  see  that  the  frequency  spectrum  of  the  estimated  mean  kernel 
has  the  salient  features  contained  in  the  experimental  pulse  echo  wavelet  measured  in 
experiment  (B).  The  center  frequency  of  the  spectra  fall  very  close  to  the  nominal  center 
frequency  of  the  transducer,  3.5  MHz.  The  measured  wavelet  has  a  slightly  narrower 
main  lobe,  probably  due  to  the  ringing  introduced  by  the  wire  target  used  in  experiment 
(B).  However,  it  should  be  kept  in  mind  that  the  measured  kernel  was  obtained  from 
simplified  underwater  experiments,  which  do  not  truly  represent  the  situation  inside  the 
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tissue  mimicking  phantom. 

The  differences  between  the  average  estimated  kernel  and  the  measured  kernel  can  be 
attributed  to: 

(El)  the  variance  associated  with  the  estimation  of  third  order  cumulants, 

(E2)  the  much  higher  dispersive  attenuation  and  aberrations  encountered  inside  the  phan¬ 
tom  in  experiment  (A),  compared  with  those  in  experiment  (B)^ 

(E3)  the  non-whiteness  of  the  underlying  scatter  response  of  the  phantom  that  violates 
assumption  (A2), 

(E4)  approximate  realization  of  a  point  target  by  a  line  target,  in  experiment  (B). 

Averaging  techniques  that  we  used  with  both  the  cumulant  estimates  of  observations  and 
estimated  kernels  themselves,  make  the  contributions  to  total  error  from  (El)  small.  The 
errors  introduced  by  approximating  a  point  target  by  a  thin  wire  target  does  not  introduce 
serious  errors  either;  this  method  is  being  used  in  the  testing  of  ultrasound  transducers  in 
research  environments,  [32],  [21]. 

The  major  contribution  to  the  differences  in  estimated  and  measured  kernels  are  from 
{E2)  and  (E3).  In  fact,  by  virtue  of  modifying  the  frequency  spectrum  of  the  received 
signal,  attenuations  appear  as  one  component  contributing  to  the  non-whiteness  (color)  of 
the  scatterer  response,  and  thereby  affects  the  kernel  estimation  procedure.  To  minimize 
the  effect  of  attenuations  on  the  stationarity  of  our  observations,  we  considered  only  short 
segments  (eg.  64  or  128  samples  at  a  13.3MHz  sampling  rate)  from  each  A-line  of  inter¬ 
est.  Identifying  and  compensating  for  the  color  of  the  tissue  response  can  be  achieved  by 
applying  the  blind  deconvolution  procedure  proposed  in  [31]  on  different  image  lines  [2]. 

The  results  we  obtained  support  our  model  assumptions  (Al)  — (A3).  It  should  be  noted 
that  ( A2)  actually  contradicts  the  Gaussianity  assumption,  which  is  commonly  made  for  the 
tissue  response.  However,  non  Gaussian  models  for  the  backscattered  ultrasonic  signals  have 
been  suggested  and  used  in  the  past  [42],  [34],  [29],  [30].  In  the  method  proposed  here,  if  the 
tissue  response  were  Gaussian,  it  would  have  been  suppressed  in  the  bispectrum  domain, 
thereby  rendering  the  estimation  of  an  axial  blurring  kernel  an  impossibility.  Hence,  our 
success  in  the  estimation  supports  the  non  Gaussianity  assumption  for  RF  ultrasound 
data.  The  structure  of  the  measured  kernel  {experiment  (BJ)  suggest  that  the  axial  kernel 
is  indeed  a  non-minimum  phase  signal;  the  method  proposed  in  this  paper  is  capable  of 
estimating  non-minimum  phase  signals  in  contrast  to  existing  techniques  which  can  only 
estimate  the  minimum-phase  equivalent  of  the  true  non-minimum  phase  signal. 
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We  consider  the  agreement  between  the  estimated  and  measured  kernels  satisfactory 
within  the  ability  of  underwater  measurements  (Measurement  B)  to  match  the  situation 
inside  the  phantom. 

6.3.2  HOS-based  distortion  estimation  and  the  deconvolution  of  B-mode  im¬ 
ages  assuming  white  tissue  response  (phantom  data,  human  liver,  human 
breast) 

A.  Images  of  a  tissue  mimicking  phantom  obtained  with  a  single  element  transducer. 
Having  established  the  validity  of  our  assumptions  and  the  kernel  estimation  technique 
in  the  axial  direction,  we  proceeded  to  estimate  lateral  distortion  kernels  from  the  image 
obtained  in  experiment  (A).,  and  to  perform  deconvolution.  The  lateral  distortion  kernels 
we  estimated  at  locations  covering  a  range  of  depths  6.58cm-9.42cm  from  the  transducer 
surface  can  be  seen  in  Fig.  3(a);  their  spectra  are  shown  in  Fig.  3(b).  In  these  estimations, 
we  used  N  =  32  and  M  =  10.  Five  adjacent  lateral  lines  were  used  in  each  kernel  estimation. 

In  the  lateral  direction,  the  estimated  kernels  have  characteristics  similar  to  the  lateral 
beam  profiles  of  the  imaging  system.  Unfortunately,  there  is  no  direct  way  to  verify  the 
results  of  our  lateral  estimations.  In  conventional  transverse  beam  profile  estimations,  the 
experiments  are  usually  done  under  water,  and  the  pressure  profiles  are  peak  detected.  This 
process  masks  instantaneous  features  of  the  beam,  hence  the  results  of  such  experiments 
can  not  be  used  in  a  verification.  Moreover,  it  has  been  shown  that  in  the  presence  of 
aberration,  the  lateral  point  spread  function  may  undergo  significant  modifications  [6], 
[39],  possibly  serious  enough  to  change  diagnoses  in  clinical  situations  [39].  Time  histories 
of  two  dimensional  instantaneous  beam  profiles  at  the  focus,  shown  in  [22],  are  in  agreement 
with  the  observation  that  the  lateral  beam  profile  can  depart  heavily  from  the  ideal  in  the 
presence  of  tissue  inhomogeneities.  These  results  suggest  that,  in  principle,  underwater 
experiments  or  theoretical  formulae  that  do  not  take  aberration  in  to  account,  can  not  be 
used  in  verifying  lateral  blurring  kernels  estimated  from  complex  targets.  Further,  it  should 
be  noted  that  even  though  our  lateral  kernels  are  estimated  at  individual  lateral  image  lines 
(i.e.,  at  a  fixed  time)  nearby  lines  contribute  to  the  data  at  a  given  time  by  virtue  of  the 
spatio-temporal  nature  of  Ultrasonic  System  Impulse  Response  defined  earlier. 

Once  the  distortion  kernels  have  been  identified,  retrieving  the  corresponding  true  tis¬ 
sue  response  becomes  a  typical  deconvolution  process.  For  the  deconvolution,  we  used  the 
constrained  Wiener  filter  technique  described  in  [40].  In  performing  the  axial  deconvolu¬ 
tion,  we  used  the  mean  HOS-based  axial  distortion  kernel,  Vm{t)  shown  in  Fig.  2.,  which 
was  obtained  with  the  method  of  Section  6.1.1;  in  the  lateral  direction,  the  mean  lateral 
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distortion  kernel  obtained  as  the  average  of  the  estimated  kernels  shown  in  Fig.  3  was  used. 

Fig.  1(b)  shows  the  logarithmically  compressed  envelope  of  the  B-scan  image,  derived 
from  the  laterally  deconvolved  RF  data  corresponding  to  the  original  image  shown  in  Fig. 
1(a).  The  axially  deconvolved  image  is  in  Fig.  1(c)  and  the  laterally  followed  by  axially 
deconvolved  image  is  in  Fig.  1(d).  Each  image  represents  an  area  of  2.5cmX2.0cm  in  the 
true  tissue  space.  For  the  sake  of  computational  simplicity,  we  assumed  that  the  lateral 
(and  axial)  blurring  kernel  do  not  vary  significantly  over  different  RF-image  lines  or  different 
image  depths  in  the  amounts  we  are  concerned  with  here.  This  allowed  us  to  use  a  single 
blurring  kernel,  in  each  of  the  axial  and  lateral  directions,  to  deconvolve  the  entire  image. 
Using  the  average  lateral  kernel  to  laterally  deconvolve  a  region  of  the  image  is  justifiable 
by  the  fact  that  all  the  kernels  estimated  from  the  depth  range  covering  the  image  shown 
here  (6.58cm  to  9.42cm)  show  a  certain  degree  of  similarity.  This  may  be  due  to  the  fact 
that  the  transducer  we  used  had  a  focal  zone  of  6- 13cm,  into  which  range  the  lateral  image 
lines  considered  above  belong.  Moreover,  the  phantom  we  used  in  experiments  had  not 
been  designed  to  simulate  strong  aberrating  effects.  As  for  the  axial  deconvolutions,  the 
pulse  echo  wavelet  (the  major  contributor  to  the  axial  distortion  kernel)  of  the  imaging 
system  is  known  to  stay  fairly  constant  over  different  RF-lines  across  the  image  [18],  [16]; 
our  results  shown  in  Fig.  2  also  support  that  hypothesis.  However,  the  method  proposed 
here  is  still  valid,  if  one  needs  to  estimate  kernels  for  localized  regions  or  individual  image 
lines.  To  capture  fine  details  in  the  pass  band  of  the  spectrum,  one  may  still  require  a  line 
by  line,  or  region  by  region,  lateral  deconvolution. 

Quantifying  resolution  improvement 

According  to  Fig.  1,  the  deconvolution  in  the  RF-domain  results  in  a  significant  reduc¬ 
tion  in  the  size  of  speckles,  suggesting  a  gain  in  resolution.  In  order  to  quantify  the  apparent 
increase  in  resolution,  we  defined  a  measure  of  resolution  based  on  the  2-D  auto-covariance 
function  of  the  image.  Auto- covariance  function  of  an  image  was  computed  on  the  RF  data 
corresponding  to  the  image,  with  the  peak  value  of  autocovariance  normalized  to  1.0.  The 
lateral  slice,  Too,  through  the  peak  of  the  2-D  auto-covariance  function  was  considered  in 
defining  the  lateral  resolution.  The  axial  slice,  Aoo,  through  the  peak  in  defining  the  axial 
resolution.  Similarly,  the  lateral  resolution,  was  defined  as  the  reciprocal  of  the  width 
of  Too,  at  d  dB  below  the  peak  of  the  slice.  Axial  resolution  was  defined  similarly  except 
for  the  fact  that  the  envelope  of  the  absolute  value  of  the  slice  Aoo  was  used.  In  this  paper 
we  used  d  —  5.00  and  d  —  10.00.  Defining  the  resolution  at  two  different  dB  levels  allows 
us  to  get  a  better  idea  of  the  shape  of  the  auto-covariance  function. 

The  reason  for  defining  the  axial  resolution  based  on  the  envelope  of  absolute  value  of 
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Aoo  is  as  follows.  We  define  the  resolution  in  terms  of  the  width,  /q,  of  the  main  lobe  of  the 
auto-covariance  function  because,  lo  is  a  measure  of  the  average  “smallness”  of  the  basic 
building  elements  of  the  image.  In  the  axial  direction,  the  representation  of  a  point  target  is 
an  oscillatory,  time-limited  signal.  Therefore,  even  when  the  target  is  a  point  in  space,  the 
received  signal  would  have  an  auto-covariance  function  which  shows  oscillatory  behavior 
(eg.  Fig.  4(d)).  In  this  case,  the  central  lobe  conveying  information  on  the  average  image 
element  size  is  the  envelope  of  the  absolute  value  of  the  auto-covariance  function. 

The  resolution  was  defined  in  the  RF  domain,  because  our  kernel  estimation  and  decon¬ 
volution  procedure  was  done  entirely  in  the  RF-domain.  The  process  of  envelope  detecting 
itself  introduces  blurring  leading  to  an  inaccurate  estimate  of  the  potential  resolution  deliv¬ 
ered  by  deconvolution,  thus  making  the  envelope  detected  signal  domain  a  less  than  ideal 
place  to  define  resolution.  Fig.  4(a)  shows  a  2D  shaded  mesh  plot  of  the  auto-covariance 
function  of  the  RF-data  corresponding  to  original  image  in  Fig.  1(a).  Fig.  4(b)  shows 
the  auto-variance  function  after  lateral  and  axial  deconvolutions.  Absolute  values  of  the 
auto-covariance  function  have  been  plotted  for  easy  visualization.  The  slice  Too  is  shown  in 
Fig.  4(c),  where  the  solid  line  indicates  the  slice  corresponding  to  the  original  image  (Fig. 
1(a))  and  the  dotted  line  that  of  the  laterally  and  axially  deconvolved  image  (Fig.  1(c)). 
Corresponding  plots  of  the  axial  slice  Aoo  are  in  Fig.  4(d). 

The  general  decrease,  in  all  directions,  of  the  width  of  the  main- lobe  of  the  auto- 
covariance  function  due  to  deconvolution  is  obvious  from  Figs.  4(a)  and  4(b).  From  Fig. 
4(c),  we  computed  the  lateral  resolution  gain,  defined  as: 
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where  1^^°^  and  respectively  represent  lateral  resolutions  before  and  after  deconvolu¬ 
tion.  A  similar  definition  holds  for  the  axial  resolution  gain.  The  lateral  resolution  gain  at 
d  =  bdB  and  d  =  lOdB  levels,  between  the  original  (Fig.  4(a))  and  the  laterally  followed 
by  axially  deconvolved  image  (Fig.  4(c))  was  found  to  be  2.7  and  3.0  respectively.  The 
corresponding  figures  for  axial  resolution  gain,  was  1.73  and  1.72  respectively. 

Our  results  indicate  that  the  lateral  resolution  gain  is  higher  than  the  axial  resolu¬ 
tion  gain.  Coupled  with  the  fact  that,  in  the  original  image  blurring  was  much  higher  in 
the  lateral  direction,  we  can  conclude  lateral  deconvolution  is  mostly  responsible  for  the 
improvement  in  overall  resolution. 

To  investigate  the  effects  of  deconvolution  on  “speckle  noise”  levels,  we  defined  the 
signal-to- noise  ratio  as,  SNR  =  (^),  where  g,  is  the  mean  of  the  image  and  cr^  is  the  variance. 
Since  we  performed  our  deconvolutions  in  the  RF  domain,  our  SNR  calculations  were  also 
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done  in  the  RF-domain,  on  absolute  values  of  RF-data.  Computed  over  the  original  image 
(Fig.  1(a)),  SNR  —  1.24;  computed  over  the  laterally  and  axially  deconvolved  image, 
SNR  =  1.26.  After  only  an  axial  deconvolution,  (Fig.  1(c)),  SNR  —  1.25,  and  after  only  a 
lateral  deconvolution  (Fig.  1(b)),  SNR  =  1.23.  Based  on  these  numbers,  we  conclude  that 
the  deconvolution  results  in  a  gain  in  resolution,  but  does  not  significantly  alter  speckle 
noise  levels.  The  results  of  SNR  computations  and  resolution  gains  have  been  summarized 
in  Tables  1  and  2. 

B.  Images  of  a  tissue  mimicking  •phantom  obtained  with  clinical  equipment 


In  order  to  investigate  the  performance  of  our  technique  with  B-scan  images  taken  from 
modern  clinical  equipment,  we  estimated  axial  and  lateral  distortion  kernels  from  the  RF- 
data  collected  in  experiment  (C).  Fig.  5(a)  shows  the  logarithmically  compressed  envelope 
of  the  original  B-scan  image.  Average  axial  and  lateral  kernels  estimated  from  the  RF-data 
corresponding  to  Fig.  5(a)  are  illustrated  in  Fig.  6(a)  and  6(b)  respectively.  These  kernels 
indicate  the  average  of  20  kernels  estimated  from  the  tissue  mimicking  background  region 
of  the  phantom,  between  the  two  target  cylinders  (see  Fig.  5).  The  spectra  of  the  average 
axial  kernel  is  shown  in  Fig.  6(c)  and  that  of  the  average  lateral  kernel  in  Fig.  6(d). 

The  result  of  lateral  deconvolution  is  shown  in  Fig.  5(b);  the  result  of  lateral  followed 
by  axial  deconvolution  is  shown  in  Fig.  5(c).  Average  estimated  kernels  shown  in  Fig.  6 
were  used  in  both  axial  and  lateral  deconvolutions.  Clearly,  the  deconvolution  has  resulted 
in  a  reduced  speckle  size  and  cleaner,  better  defined  boundaries  of  the  target  cylinders.  It 
is  also  evident  that  the  attenuations  associated  with  the  imaging  process  show  up  much 
clearly  in  the  deconvolved  image.  This  may  have  significance  in  clinical  imaging  situations, 
where  attenuation  properties  of  tissue  convey  important  diagnostic  information. 

Fig.  7(a)  and  7(b)  illustrate  the  absolute  value  of  the  auto-covariance  function  of  the 
RF-data  corresponding  to  the  original  (Fig.  7(a))  image  and  the  laterally  and  axially 
deconvolved  image  (Fig.  7(d)),  respectively.  Fig.  7(c)  and  7(d)  show  the  lateral  and 
axial  slices  Loo  and  Aoo  used  to  compute  resolution  gains.  Based  on  data  corresponding  to 
Figs.  5  and  7,  speckle  SNR  and  resolution  gains  were  computed  and  tabulated  in  Tables 
1  and  2.  Figs.  5,  7  and  Tables  1  and  2  lead  to  the  conclusion  that  considerable  resolution 
enhancement  is  possible  with  deconvolution,  and  that  the  process  of  deconvolution  does  not 
affect  the  speckle  noise  significantly.  Once  again,  lateral  resolution  is  found  to  be  mostly 
responsible  for  the  overall  improvement  in  the  image. 

C.  Clinical  images  of  a  human  liver 

To  evaluate  our  method  with  clinical  images,  we  obtained  a  B-scan  of  a  liver  image  as 
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described  in  experiment  (D),  section  3.1.  Figure  8(a)  shows  a  part  of  the  logarithmically 
compressed  envelope  of  the  original  image.  Fig  9(a)  and  9(b)  show  the  average  axial  and 
lateral  distortions  estimated  from  the  RF-data  corresponding  to  Fig.  8(a).  Their  spectra 
can  be  seen  in  Figs.  9(c)  and  9(d). 

The  results  of  lateral  deconvolution  is  shown  in  Fig.  8(b);  results  of  lateral  followed  by 
axial  is  in  Fig.  8(c).  In  all  cases,  logarithmically  compressed  envelope  of  the  images  has 
been  displayed. 

Again,  to  visualize  the  improvement  in  resolution,  we  plotted  the  auto-covariances  of 
the  original  and  deconvolved  images.  Fig.  10(a)  and  10(b)  illustrates  the  shaded  2-D  mesh 
plot  of  the  auto-covariances  of  RF-data  corresponding  to  Figs.  8(a)  and  8(c)  respectively. 
Corresponding  Lqo  and  Aqq  slices  are  shown  in  Figs.  10(c)  and  10(d).  Slices  corresponding 
to  original  image  are  shown  in  solid  lines,  whereas  those  of  the  deconvolved  image  are  shown 
in  dotted  lines. 

The  lateral  resolution  gain,  at  d  =  5dB  and  d  —  lOdR  was  found  to  be  3.0  and  2.5 
respectively.  Corresponding  numbers  for  the  axial  resolution  gain  were  1.7  and  1.5.  Speckle 
SNR  ratio  computed  for  the  original  image  (Fig.  8(a))  was  SNR  =  1.20;  after  lateral  and 
axial  deconvolutions  (Fig.  8(c))  SNR  =  1.33.  Tables  1  and  2  summarize  these  results. 

These  results  suggest  that  considerable  gain  in  resolution  is  possible  with  axial  and 
lateral  deconvolution  of  clinical  images. 

Clinical  ultrasound  images  of  the  breast 

Results  obtained  in  experiment  E  described  in  Section  6.2  are  shown  here. 

Fig.  11(a)  shows  the  B-scan  breast  image  of  Patient  1.  The  result  of  lateral  followed 
by  axial  image  deconvolution  is  shown  in  Fig.  11(b).  In  both  cases,  the  logarithmically 
compressed  envelope  of  the  RF-data  has  been  used  for  display  purposes.  Each  image 
represents  an  axial  depth  of  about  2.5  cm  in  the  true  tissue  space. 

According  to  Fig.  11,  the  deconvolution  in  the  RF-domain  has  resulted  in  a  significant 
reduction  in  the  size  of  speckles,  suggesting  a  gain  in  resolution.  In  order  to  quantify  the 
apparent  increase  in  resolution  we  computed  the  resolution  gain  for  the  image,  before  and 
after  deconvolution.  At  the  10  dB  and  5  dB  levels,  the  resolution  gains  were  found  to  be 
2.0  and  2.6  (see  Table  3).  The  deconvolution  process  lead  an  improvement  in  the  image 
speckle  SNR  ratio,  from  SNR  =  0.68  to  0.997  (see  Table  3).  These  objective  measures 
suggest  that  the  deconvolution  deconvolution  process  is  capable  of  improving  the  resolution 
of  breast  ultrasound  images  without  compromising  the  speckle  SNR. 

Similar  results  for  the  breast  images  of  Patients  2-5  are  shown  in  Figs.  11  to  15  and 
Tables  3,4.  All  these  results  consistently  point  to  the  fact  that  the  deconvolution  leads  to 
an  improvement  in  the  resolution. 
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However,  the  real  improvement/deterioration  due  to  deconvolution  can  be  fully  appre¬ 
ciated  only  by  trained  end  users,  i.e.  radiologists,  who  make  diagnostic  decisions  based 
on  the  images.  As  a  preliminary  means  of  obtaining  the  radiologists  response,  we  showed 
the  original  and  deconvolved  images  to  a  radiologist  at  the  Breast  Imaging  Center  of  the 
Thomas  Jefferson  University  Hospital,  Philadelphia.  The  responses  are  as  follows: 

•  The  radiologist  consistently  preferred  the  deconvolved  images  over  the  unprocessed 
ones  in  all  the  cases,  on  the  grounds  that  finer  details  can  be  easier  to  perceive  from 
the  deconvolved  image  due  to  its  smaller  speckle  sizes  (fine-grainess). 

•  In  the  case  of  Patient  i,  (see  Fig.  11),  the  tumor  (indicated  by  arrows  in  Fig.  11(a)) 
and  its  spread  (boundaries)  were  easier  to  detect/determine  from  the  deconvolved 
image. 

•  In  the  case  of  Patient  3,  (see  Fig.  13),  the  deconvolution  resulted  in  a  better  visu¬ 
alization  of  the  jaggerdness  of  the  tumor  boundary,  disclosing  tumor  infiltration  to 
surrounding  tissue  and  probably  suggesting  the  malignant  nature  of  the  tumor.  This 
detail  was  not  easily  perceived  in  the  unprocessed  image. 

•  In  the  case  of  Patient  4  ,  (see  Fig.  14),  the  smooth  nature  of  the  boundary  of  the 
tumor  has  become  very  easy  to  perceive,  compared  with  the  unprocessed  image. 
Smooth  boundaries  may  suggest  a  benign  tumor. 

•  In  the  case  of  Patient  5  ,  (see  Fig.  15),  the  finer  details  have  become  easier  to  perceive 
in  the  deconvolved  image.  More  importantly,  deconvolution  resulted  in  one  region 
(indicated  by  arrows)  becoming  a  suspected  area.  It  was  not  very  easy  to  judge  in 
the  first  image. 

•  Deconvolution,  even  though  makes  the  image  texture  finer,  does  not  lead  to  a  dramatic 
change  in  the  nature  of  clinical  ultrasound  images  familiar  to  radiologists.  As  a  result, 
it  is  not  hard  to  get  used  to  deconvolved  images  in  diagnosis. 

•  It  would  be  very  useful  if  the  deconvolution  can  be  done  real-time.  One  of  the  major 
advantages  of  medical  ultrasound  is  its  real-time  in  nature,  which  results  in  dynamic 
images  conveying  more  information  than  a  still  picture  scan. 

Although  these  responses  are  preliminary,  they  are  very  positive  and  encourages  one  to 
conduct  a  systematic  analysis  into  the  advantages  of  image  deconvolution  in  clinical  diag¬ 
nosis.  This  can  be  attempted  ones  a  large  number  of  images  are  available  for  comparison. 
A  large  number  of  radiologists  should  also  be  involved  in  the  process. 
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6.3.3  HOS-based  estimated  distortions  and  deconvolution  of  B-mode  images 
assuming  non- white  tissue  response 

To  investigate  the  significance/implications  of  the  whiteness  assumption  on  the  tissue  re¬ 
sponse  we  attempted  estimating  lateral  and  axial  kernels  from  a  clinical  B-scan  image 
described  in  Section  6.2,  Experiment  D.  A  part  of  the  image  approximately  7cmX2cm  in 
size  (in  true  tissue  space)  and  containing  a  hypoechoic  metastasis  of  the  liver,  is  shown  in 
Fig.  16(a).  Logarithm  of  the  envelope  has  been  used  for  display  purposes. 

Axial  Distortion  Kernels 

Some  typical  examples  of  the  combined  axial  kernels  fc“(n)  estimated  from  the  corre¬ 
sponding  rf-image  are  shown  in  Fig.  17.  In  order  to  reduce  the  estimation  variance  we 
concatenated  3  adjacent  rf  A-lines  in  each  case.  The  effect  of  non-stationarities  was  min¬ 
imized  by  considering  short  segments  of  data,  128  samples  in  this  case,  from  each  A-line. 
The  kernels  shown  here  correspond  to  a  region  of  approximately  12mm^,  in  the  original 
tissue  space.  Figs.  18(a)  and  (b)  show  two  examples  of  the  kernels  estimated  from  within 
the  tumor,  i.e.  kii{n)^  and  Fig.  18(c)  and  (d)  show  two  kernels  estimated  from  outside  the 
tumor,  i.e.  kio{n).  The  corresponding  spectra  are  illustrated  in  Fig.  19;  the  spectra  of  the 
rf  A-lines  used  in  estimating  the  kernels  and  those  of  the  kernels  themselves  are  respectively 
shown  by  dotted  and  solid  lines. 

Several  observations  are  in  order  now:  (a)  Fig.  18  depicts  that  there  is  a  good  match 
between  the  spectra  of  combined  kernels  and  corresponding  rf  data.  Figs.  17(a)  and(b) 
show  that  there  is  a  similarity  between  combined  kernels  k^Q{n)  estimated  outside  the 
tumor;  the  same  is  true  for  combined  kernels  kfj[{n)  estimated  inside  the  tumor,  and  (c) 
kfoin)  are  slightly  different  from 

The  difference  between  kfj{n)  and  kfQ{n)  can  be  attributed  to  the  differences  between 
color  of  the  tissue  responses  occurring  inside  and  outside  the  tumor.  Since  tumors  change 
the  structure  of  normal  tissues  by  definition,  indeed  we  expect  to  see  a  change  in  the  color 
of  tissue  response. 

Applying  the  blind  deconvolution  scheme  developed  in  [31]  on  the  two  combined  axial 
kernels  k^j{n)  and  klQ{n)  which  are  shown  in  Figs.  17(a)  and  (c)  respectively,  we  identified 
the  color  of  the  tissue  responses  inside  the  tumor  (tij(n))  and  outside  the  tumor  (^io(^))- 
Obtained  color  terms  tJo(n)  and  tj/(n)  are  shown  in  Fig.  19,  together  with  their  spectra. 
It  is  evident  that  the  color  of  the  tissue  response  outside  the  tumor  (Fig.  19(b), (d))  is 
different  from  the  one  inside  the  tumor  (Fig.  19(a)  and  (c)).  Apparently,  this  reflects  the 
change  in  the  tissue  structure  inside  the  tumor. 

Once  we  know  the  color  of  the  tissue  response  t“(n),  we  can  estimate  the  true  axial 
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distortion  kernel  ft“(n)  from  eq.(lO)  and  (11),  through  a  deconvolution  process.  We  used 
the  constrained  Wiener  filter  technique  described  in  [40]  for  this  purpose.  Fig.  20(a) 
illustrates  the  resulting  estima,tes  for  h°'{n).  The  corresponding  spectra  are  plotted  in  Fig. 
20(b). 

Lateral  Distortion  Kernels 

In  Fig.  21,  two  of  the  typical  combined  lateral  distortion  kernels  estimated  inside  and 
outside  of  the  tumor  are  shown,  together  with  their  spectra.  In  Figs.  21(b)  and  (c), 
solid  lines  indicate  the  spectra  of  data  which  were  used  to  estimate  the  distortion  kernels, 
and  dotted  lines  indicate  the  spectra  of  the  kernels  themselves.  In  order  to  account  for 
depth  related  beam  spreading,  both  kernels  were  estimated  at  the  same  axial  depth,  5.0cm, 
in  the  true  tissue  space.  As  before,  20  adjacent  lateral  data  lines  in  the  rf-image  were 
concatenated  to  form  the  basis  for  the  estimations.  Each  lateral  kernel  represents  an  area 
of  17mm^  (approx.)  in  the  original  tissue  space.  From  Fig.  21,  we  can  see  that  the  lateral 
distortion  kernels  estimated  from  inside  and  outside  of  the  tumor  assume  two  slightly 
different  shapes.  The  same  observation  is  true  for  the  spectra  of  the  rf-data  lines  used  to 
estimate  these  kernels. 

Using  the  same  procedure  followed  in  the  case  of  axial  distortion  kernels,  we  estimated 
the  color  of  the  tissue  response  (lateral  direction),  working  on  the  lateral  distortion  kernels 
shown  in  Fig.  21.  The  resulting  terms  corresponding  to  the  color  of  the  tissue  response 
are  shown  in  Fig.  22,  together  with  their  spectra.  It  is  clear  that  the  color  of  the  tissue 
response,  t\{n)  is  different  in  the  two  regions  of  the  image  considered. 

Following  the  same  procedure  as  in  the  case  of  axial  distortion  kernels,  we  estimated 
the  true  lateral  distortion  kernel.  Resulting  color-free  lateral  distortion  kernels,  h\n)  and 
their  spectra  are  shown  in  Fig.  23. 

Using  the  lateral  distortion  kernel  h\n)  estimated  above,  we  deconvolved  the  rf-data 
corresponding  to  the  image  shown  in  Fig.  16(a).  The  result  of  lateral  deconvolution  is 
shown  in  Fig.  16(b),  where  the  logarithmically  compressed  envelope  has  been  used  for 
display  purposes.  The  result  of  lateral  deconvolution  using  the  kernel  h\n)^  followed  by 
the  axial  deconvolution  using  h°^(n)  estimated  earlier,  is  in  Fig.  16(c).  Note  that  we  used 
only  a  single  lateral  kernel  to  deconvolve  the  entire  image.  This  may  be  justifiable  as  a  first 
approximation,  considering  the  depth  (approx.  2  cm  in  the  real  tissue  space)  of  the  image 
involved. 

Deconvolution  of  the  image  was  also  performed  using  the  combined  distortion  kernels. 
These  kernels  are  basically  the  same  kernels  that  would  result  under  the  white  assumption 
about  the  tissue  response.  The  deconvolution  results  were  identical  to  the  ones  shown  in 
Figs.  16(b)  and  (c). 
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Based  on  Fig.  16,  we  can  conclude  that  there  is  a  considerable  resolution  enhancement, 
especially  after  the  lateral  deconvolution,  while  the  resolution  improvement  is  the  same 
for  both  combined  and  true  distortion  kernels.  The  latter  observation  leads  us  to  the 
conclusion  that  tissue  response  can  be  safely  assumed  white  for  deconvolution  purposes. 
Although  there  is  some  difference  in  the  combined  kernels  inside  and  outside  the  tumor,  or 
otherwise,  when  tissue  structure  variations  occur,  this  difference  is  not  significant  to  affect 
the  deconvolution  result.  This  difference,  however,  emerges  as  potential  tissue  signature  [2]. 
If  we  are  interested  in  detecting  tissue  structure  changes,  (e.g.,  changes  between  tumorous 
and  tumor-free  areas)  the  combined  kernels  can  not  be  used.  The  distortion  kernels  are 
the  dominant  components  in  the  combined  kernels  and  mask  the  small  contribution  of  the 
color  of  the  tissue  response  which  carries  the  discriminating  information.  The  color  of  the 
tissue  response  as  a  tissue  signature  has  been  investigated  in  [2]  (see  also  Appendix  A). 

6.3.4  SOS-based  estimated  distortions  and  deconvolution  of  B-mode  images 

From  the  images  obtained  at  experiments  A,  C  and  D,  SOS-based  distortions  were  estimated 
with  the  method  of  Section  6.1.3.  The  images  were  subsequently  deconvolved,  and  the 
deconvolution  results  are  shown  in  Figs.  27,  28  and  29,  respectively.  The  corresponding 
resolution  gains  and  speckle  SNR  are  shown  in  Tables  5  and  6. 

Comparing  Tables  1  and  2  to  5  and  6,  we  can  conclude  that  the  resolution  improvement, 
both  axial  and  lateral,  was  consistently  superior  when  deconvolution  was  performed  with 
the  HOS-based  estimated  distortions.  Deconvolution  with  HOS-based  estimates  led  to 
axial  resolution  improvement  1.15  -  1.9  times  higher  than  deconvolution  with  SOS-based 
estimates.  The  improvement  of  lateral  resolution  when  deconvolving  using  HOS  based 
estimates  was  more  dramatic;  it  was  better  than  the  SOS  based  deconvolution  by  a  factor 
in  the  range  1.68-4.5.  The  superiority  of  the  HOS-based  deconvolution  is  also  evident 
by  comparing  the  corresponding  images  deconvolved  with  HOS  and  SOS-based  distortion 
estimates. 

7  CONCLUSIONS/FUTURE  WORK 

Processing  ultrasound  images  with  higher-order  spectral  operations  we  were  able  to  identify 
the  distortion  introduced  by  the  ultrasonic  system  and  the  medium.  The  proposed  method 
makes  it  possible  to  estimate  both  axial  and  instantaneous  lateral  blurring  kernels,  working 
on  B-mode  RF  data.  Distortion  identification  and  subsequent  cancelling  (deconvolution) 
operations  were  carried  out  on  1-D  lines  of  the  RF  image,  thereby  obviating  the  theoretical 
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difficulties  faced  by  earlier  attempts  at  beam  deconvolution  on  envelope  detected  images. 
The  method  is  capable  of  estimating  mixed  phase  distortion  kernels,  and  is  immune  to 
additive  Gaussian  noise. 

Performing  underwater  experiments  employing  single  element  transducers,  commercial 
tissue  mimicking  phantoms  and  simulated  point  targets,  we  showed  that  our  kernel  estima¬ 
tion  procedure  could  be  done  with  reasonable  accuracy.  The  accuracy  of  the  estimations 
were  verified  by  measuring  axial  kernels  underwater. 

Deconvolution  results  obtained  with  phantom  data  and  clinical  images  indicate  con¬ 
siderable  resolution  improvement.  Lateral  deconvolution  contributes  heavily  to  the  gain 
in  resolution,  where  the  resolution  was  defined  in  terms  of  the  dimensions  of  the  auto¬ 
covariance  function  of  the  image.  These  results  are  significant  because  in  unprocessed 
images,  lateral  distortions  are  known  to  be  as  much  as  3-5  times  severe  than  axial  distor¬ 
tions,  leading  to  a  change  in  appearance  in  clinical  images  depending  on  the  orientation  of 
the  scanner.  Through  lateral  deconvolution,  one  can  compensate  for  the  lateral  distortions 
and  try  to  achieve  consistent  images  independent  from  the  angular  position  of  the  scanner. 
In  the  past,  distortion  estimation  was  carried  out  exclusively  in  the  second-order  statistics 
domain.  We  demonstrated  that  deconvolution  based  on  HOS-based  distortion  estimates 
leads  to  superior  axial  and  in  particular  lateral  resolution  improvement  than  the  one  with 
SOS-based  estimates. 

Liver  tissue  is  relatively  homogeneous,  thus  the  resolution  improvement  via  deconvo¬ 
lution  came  at  little  surprise.  The  highly  inhomogeneous  nature  of  breast  tissue,  on  the 
other  hand,  made  the  application  of  the  proposed  methodology  to  breast  data  a  not  at  all 
straightforward  task.  It  was  really  exciting  to  see  that  indeed  the  proposed  methodology 
led  to  significant  resolution  gains  when  applied  on  human  breast  images.  The  resolution 
gains  were  quantified  based  on  our  standard  mathematical  measures,  but  most  important 
based  on  the  opinion  of  our  medical  collaborators  at  Thomas  Jefferson  University  Hospital. 
Trained  radiologists  compared  the  images  before  and  after  processing  and  found  that  the 
processed  images  revealed  fine  details  that  were  hidden  in  the  original  images,  and  which 
could  be  critical  in  determining  the  malignant  or  not  nature  of  the  tumor.  Although  more 
extensive  studies  need  to  be  conducted,  and  on  line  processing  should  be  considered,  we 
view  the  obtained  results  as  very  encouraging. 

It  has  been  reported  that  high  transducer  frequencies  required  for  high  resolution  imag¬ 
ing  actually  lead  to  lower  resolution  in  the  presence  of  increased  aberration  in  breast  at 
higher  frequencies  [8].  In  [36]  it  is  pointed  out  that  medium  inhomogeneities  are  also  impor¬ 
tant  in  systems  where  a  larger  aperture  is  used  for  higher  lateral  resolution.  In  general,  the 
aberration  correction  is  of  fundamental  importance  in  high  resolution  medical  ultrasonic 
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systems.  The  ability  of  the  proposed  method  to  compensate  for  convolution  components 
of  aberration  holds  promise  in  improving  the  diagnostic  value  of  B-mode  breast  images 
beyond  the  capabilities  of  hardware. 

Deconvolution,  viewed  as  an  image  de-blurring/restoration  operation,  should  reveal 
those  small  structures  (such  as  early  stage  tumors)  that  had  been  hidden  away  from 
view  due  to  imaging  distortions  (including  convolutional  components  of  aberrations),  but 
otherwise  would  show  up  on  a  hypothetical,  perfect  imager.  A  deconvolved  image  also 
means  access  to  a  distortion-free  tissue  signal,  which  is  largely  independent  of  the  imaging 
system.  Therefore,  tissue  characterization  schemes  which  ideally  require  imaging-system- 
independent  data  could  be  based  on  deconvolved  RF  images.  In  addition,  the  estimated 
distortion  kernels  themselves  carry  information  on  the  statistical  structure  of  scatterer  field 
as  manifested  through  the  color  of  the  scatterer  response  (tissue  response),  attenuations 
and,  propagation  non-linearities  associated  with  the  phantom  (tissue).  Indeed  our  prelimi¬ 
nary  investigation  suggested  that  the  color  of  the  tissue  response  is  an  excellent  candidate 
as  a  tissue  signature. 
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9  FIGURES /TABLES 


Table  1:  Resolution  gains  due  to  deconvolution  with  HOS-based  estimated  distrotions 


Resolution  Gains 

lateral 

axial 

5dB 

lOdB 

5dB 

lOdB 

Experiment  (A) 

2.7 

3.0 

n 

1.7 

Experiment  (C) 

5.2 

4.2 

1.9 

Experiment  (D) 

2.5 

1.7 

1.5 

Table  2:  The  effect  of  deconvolution  with  HOS-based  distortions  on  the  speckle  SNR 
ratio  _ _ _ _ _ 


SNR 

original 

deconvolved 

Experiment  (A) 

1.24 

1.26 

Experiment  (C) 

1.17 

1.10 

Experiment  (D) 

1.20 

1.30 

Table  3:  Resolution  gains  in  breast  images  due  to  HOS-based  deconvo¬ 
lution.  _ _ 


Resolution  Gains 

lateral 

axial 

5dB 

lOdB 

5dB 

lOdB 

Patient  1 

5.6 

2.0 

8.6 

5.6 

Patient  2 

4.8 

4.0 

1.9 

1.8 

Patient  3 

3.0 

2.1 

1.7 

1.6 

Patient  4 

2.2 

1.8 

2.1 

1.9 

Patient  5 

2.9 

1.9 

1.3 

1.2 
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Table  4:  The  effect  of  HOS-based  deconvolution  of  breast  images  on  the 
speckle  SNR  ratio. _ 


SNR 

original 

deconvolved 

Patient  1 

0.68 

0.99 

Patient  2 

0.46 

0.59 

Patient  3 

0.74 

0.84 

Patient  4 

0.51 

0.58 

Patient  5 

0.75 

0.84 

Table  5:  The  elFect  of  deconvolution  with  SOS-based  distortions  on  the  speckle  SNR 
ratio 


Resolution  Gains 

lateral 

axial 

5dB 

lOdB 

5dB 

lOdB 

Experiment  (A) 

1.6 

1.35 

1.53 

1.20 

Experiment  (C) 

1.85 

0.93 

1.48 

1.64 

Experiment  (D) 

1.50 

1.44 

1.50 

1.23 

Table  6:  The  effect  of  deconvolution  with  SOS-based  distortions  on  the  speckle  SNR 
ratio  _ 


SNR 

original 

deconvolved 

Experiment  (A) 

1.21 

1.24 

Experiment  (C) 

1.2 

1.20 

Experiment  (D) 

1.10 

1.17 
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(c) 


(d) 


Figure  1:  A  speckle-only  part  of  the  ultrasound  image  of  the  tissue  mimicking  phantom, 
obtained  with  a  focused  single  element  transducer  experiment  (A);  (b)  the  result  of  late 
ral  deconvolution;  (c)  the  result  of  axial  deconvolution  and,  (d)  the  result  of  lateral  fol¬ 
lowed  by  axial  deconvol  ution.  In  all  cases,  the  logarithmically  compressed  envelope  is 
shown.  Deconvolution  was  performed  with  the  higher-order  statistics  based  estimated 
distortions. 
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Frequency  /  Hz  x  10^ 


Figure  2:  Underwater  experiments  {experiments  (A)  and  (B)):  (a)  axial  kernels  estimated 
at  various  lateral  positions  in  the  B-mode  image  of  the  tissue  mimicking  phantom  (dotted 
lines);  average  of  the  estimated  kernels  (solid  line).  The  experimental  kernel,  measured 
as  the  reflection  off  a  0.812mm  diameter  wire  surface,  under  water.  The  close  agreement 
between  the  estimated  and  the  measured  kernels  (with  in  limits  of  experimental  errors), 
indicates  the  success  of  the  estimation  procedure,  (b)  Spectra  of  the  average  estimated 
kernel  (solid  line)  and  the  measured  kernel  (dashed  line).  The  center  frequency  of  both 
kernels  are  compatible  with  the  nominal  center  frequency  of  the  transducer,  3.5MHz. 
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Figure  3:  (a)  Lateral  distortion  kernels  estimated  at  various  axial  depths  from  the  RF  data 
corresponding  to  the  image  shown  in  Fig.  1(a).  All  the  kernels  are  estimated  within  the 
focal  zone  of  the  transducer,  6  —  13cm.  (b)  Spectra  of  the  lateral  kernels  shown  in  (a). 
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Figure  4:  Shaded  auto-covariance  function  of  the  RF-data  corresponding  to:  (a)  the  origi¬ 
nal  image  shown  in  Fig.  1(a);  (b)  the  laterally  and  axially  deconvolved  image  shown  in  Fig. 
1(d).  Absolute  values  of  the  auto-covariance  functions  have  been  shown  for  easy  visualiza¬ 
tion.  (c)  The  lateral  slice  Lqo  and  (d)  the  axial  slice  Aqq  of  the  auto-covariance  function 
corresponding  to:  the  original  image  (solid  lines),  the  deconvolved  image  (dotted  lines). 
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Figure  5:  {Experiment  (C)):  (a)  The  original  image  of  the  tissue  mimicking  phantom  ob¬ 
tained  with  the  linear  array  transducer  on  a  clinical  imaging  system;  (b)  the  result  of  lateral 
deconvolution  and  (c)  the  result  of  lateral  and  axial  deconvolution.  The  logarithmically 
compressed  envelope  has  been  used  for  display.  Deconvolution  was  performed  with  the 
HOS-based  estimated  distortions. 


Sample  Number  (lateral) 


Figure  6:  {Experiment  (C)):  (a)  The  average  of  the  estimated  axial  kernels  from  RF-data 
corresponding  to  image  shown  in  Fig.  5(a),  and  (b)  its  spectrum,  (c)  The  average  of  the 
lateral  kernels  estimated  from  RF-data  corresponding  to  Fig.  5(a)  and  (d)  its  spectrum. 
These  kernels  were  estimated  from  the  region  between  the  two  cylindrical  targets  in  Fig. 
5(a). 
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ORIGINAL 


DECONVOLVED 


LATERAL 


AXIAL 


Figure  7:  Shaded  auto-covariance  function  of  the  RF-data  corresponding  to:  (a)  the  origi¬ 
nal  image  shown  in  Fig.  5(a);  (b)  the  laterally  and  axially  deconvolved  image  shown  in  Fig. 
5(d).  Absolute  values  of  the  auto- covariance  functions  have  been  shown  for  easy  visualiza¬ 
tion.  (c)  The  lateral  slice  Too  and  (d)  the  axial  slice  Aoo  of  the  auto-covariance  function 
corresponding  to:  the  original  image  (solid  lines),  the  deconvolved  image  (dotted  lines). 
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Figure  8:  {Experiment  (D)):  (a)  A  part  of  the  clinical  image  of  a  human  liver  containing  a 
tumor;  (b)  the  result  of  lateral  deconvolution  and  (c)  the  result  of  lateral  and  axial  decon¬ 
volution.  In  all  cases  the  logarithm  of  the  envelope  is  shown.  (HOS-based  deconvolution  is 
shown  here). 
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Figure  9;  {Experiment  (D))\  (a)  The  average  of  the  estimated  axial  kernels  from  RF-data 
corresponding  to  image  shown  in  Fig.  8(a),  and  (b)  its  spectrum,  (c)  The  average  of  the 
lateral  kernels  estimated  from  RF-data  corresponding  to  Fig.  8(a)  and  (d)  its  spectrum. 
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ORIGINAL 


DECONVOLVED 


Figure  10:  Shaded  auto-covariance  function  of  the  RF-data  corresponding  to:  (a)  the 
original  image  shown  in  Fig.  8(a);  (b)  the  laterally  and  axially  deconvolved  image  shown 
in  Fig.  8(d).  Absolute  values  of  the  auto- covariance  functions  have  been  shown  for  easy 
visualization,  (c)  The  lateral  slice  Too  and  (d)  the  axial  slice  Aqo  of  the  auto-covariance 
function  corresponding  to:  the  original  image  (solid  lines),  the  deconvolved  image  (dotted 
lines). 
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Figure  11:  (Patient  1):  (a)  A  part  of  the  clinical  image  of  a  human  breast  containing  a 
tumor;  (b)  the  result  of  lateral  followed  by  axial  deconvolution.  In  all  cases  the  logarithm 
of  the  envelope  is  shown. 
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Figure  12:  {Patient  2  ):  (a)  A  part  of  the  clinical  i  mage  of  a  human  breast  containing  a 
tumor;  (b)  the  result  of  lateral  followed  by  axial  deconvolu  tion.  In  all  cases  the  logarithm 
of  the  envelope  is  shown. 


46 


(a)  (b) 


Figure  13:  {Patient  3  ):  (a)  A  part  of  the  clinical  image  of  a  human  breast  containing  a 
tumor;  (b)  the  result  of  lateral  followed  by  axial  deconvolution.  In  all  cases  the  logarithm 
of  the  envelope  is  shown. 
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(a)  (b) 


Figure  14:  (Patient  4)'  (a)  A  part  of  the  clinical  image  of  a  human  breast  containing  a 
tumor;  (b)  the  result  of  lateral  followed  by  axial  deconvolution.  In  all  cases  the  logarithm 
of  the  envelope  is  shown. 
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(a)  (b) 


Figure  15:  {Patient  5)\  (a)  A  part  of  the  clinical  image  of  a  human  breast  containing  a 
tumor;  (b)  the  result  of  lateral  followed  by  axial  deconvolution.  In  all  cases  the  logarithm 
of  the  envelope  is  shown. 
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Figure  16:  Clincal  image  of  a  human  liver,  (a)  Envelope  of  a  part  of  the  original  image 
containing  hypo-echoic  liver  metastasis  tumor;  (b)  laterally  deconvolved  image;  (c)  laterally 
followed  by  axially  deconvolved  image,  logarithm  of  the  envelope  is  shown. 
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Figure  17:  (a),(b)  Combined  axial  distortion  kernels  estimated  from  small  regions  inside 
the  liver  tumor;  (c),(d)  combined  axial  distortion  kernels  estimated  outside  the  liver  tumor. 


51 


Figure  18:  (a),(b)  Spectral  of  combined  axial  distortion  kernels  estimated  from  small  regions 
inside  the  liver  tumor;  (c),(d)  spectra  of  combined  axial  distortion  kernels  estimated  outside 
the  liver  tumor.  Solid  lines  indicate  the  spectra  of  estimated  kernels  and  dotted  lines  those 
of  the  rf  data  based  on  which  the  kernels  were  estimated. 
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Figure  19;  Color  of  the  axial  tissue  response  estimated  from  small  regions  (a)  inside  the 
tumor  and  (b)  outside  the  tumor,  (c)  Spectrum  of  the  kernel  shown  in  (a);  (d)  spectrum 
of  the  kernel  shown  in  (b). 
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Spectrum,  Unear  Scale 


Figure  20:  (a)  True  axial  distortion  kernels  reconstructed  based  on  the  combined  distortion 
kernels  estimated  from  inside  (solid  line)  and  outside  (dotted  line)  the  tumor,  and  (b)  the 
corresponding  spectra. 


Figure  21:  (a)  Combined  lateral  distortion  kernels  estimated  (a)  inside  the  tumor,  and  (b) 
its  spectrum.  Combined  lateral  distortion  kernels  estimated  (c)  outside  the  tumor,  and  (d) 
its  spectrum.  Solid  lines  indicate  the  spectra  of  estimated  kernels  and  dotted  lines  those  of 
the  rf  data  based  on  which  the  kernels  were  estimated. 
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Figure  22:  (a)  The  color  of  the  lateral  tissue  response  estimated  (a)  inside  the  tumor  and 
(b)  outside  the  tumor,  (c),  (d)  The  spectra  corresponding  to  (a)  and  (b). 
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Figure  23:  (a)  True  lateral  distortion  kernels  reconstructed  from  the  combined  kernels  esti¬ 
mated  from  inside  (dotted  line)  and  outside  (solid  line)  the  tumor,  (b)  The  corresponding 
spectra. 
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Figure  24:  A  speckle-only  part  of  the  ultrasound  image  of  the  tissue  mimicking  phan¬ 
tom,  obtained  with  a  focused  single  element  transducer  experiment  fAJ;  (b)  the  result  of 
lateral  deconvolution;  (c)  the  result  of  axial  deconvolution  and,  (d)  the  result  of  lateral 
followed  by  axial  deconvolution.  In  all  cases,  the  logarithmically  compressed  envelope  is 
shown.  Deconvolution  was  performed  with  the  second-order  statistics  based  estimated 
distortions. 
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Figure  25:  {Experiment  (C)):  (a)  The  original  image  of  the  tissue  mimicking  phantom  ob¬ 
tained  with  the  linear  array  transducer  on  a  clinical  imaging  system;  (b)  the  result  of  lateral 
deconvolution  and  (c)  the  result  of  lateral  and  axial  deconvolution.  The  logarithmically 
compressed  envelope  has  been  used  for  display.  Deconvolution  was  performed  with  the 
second-order  statistics  based  estimated  distortions. 


laterally  and  axially  deconvolved  image 


Figure  26:  {Experiment  (D)):  (a)  A  part  of  the  clinical  image  of  a  human  liver  containing 
a  tumor;  (b)  the  result  of  lateral  deconvolution  and  (c)  the  result  of  lateral  and  axial 
deconvolution.  In  all  cases  the  logarithm  of  the  envelope  is  shown,  (second-order  statistics 
based  deconvolution  is  shown  here). 
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MEETING  ABSTRACTS 


SPIE  International  Symposium  on  Optics,  Imaging  and  Instrumentation,  San  Diego,  July 
1994. 


Blind  Deconvolution  of  Ultrasound  Images 

Udantha  R.  Aheyratne,  Athina  P.  Petropulu*  and  John  M.  Reid 

Biomedical  Engineering  and  Science  Institute 
*  Department  of  Electrical  and  Computer  Engineering, 

Drexel  University,  Philadelphia,  PA  19104 

We  address  the  problem  of  improving  the  resolution  of  ultrasound  images  using  blind 
deconvolution.  The  transducer  measurement  that  forms  the  ultrasound  image  can  be  ex¬ 
pressed  as  the  convolution  of  two  terms,  the  tissue  response  and  the  ultrasonic  system 
response,  plus  additive  noise.  The  quantity  of  interest  is  the  tissue  response,  however,  due 
to  the  convolution  operation,  we  measure  a  blurred  version  of  it,  which  obscures  the  fine 
details  in  the  image.  Our  goal  is  to  remove  the  blurring  caused  by  the  ultrasonic  system, 
in  order  to  enhance  the  diagnostic  quality  of  the  ultrasound  image.  Under  the  assump¬ 
tion  that  speckle  behaves  as  an  i.i.d.  zero-mean  non-Gaussian  process,  we  were  able  to 
reconstruct  the  blurring  kernel  using  bicepstrum  operations  on  corresponding  A-scan  lines. 
Based  on  the  estimated  blurring  kernel  we  performed  deconvolution  on  the  measured  im¬ 
age.  Preliminary  results  obtained  from  ultrasound  images  of  a  tissue  mimicking  phantom 
indicate  significant  resolution  improvement. 


20^^  International  Symposium  on  Ultrasonic  Imaging  and  Tissue  Characterization,  Arling¬ 
ton,  VA,  June  1995. 

Estimating  Imaging-Distortions  and  the  Color  of  the  Tissue  Response  From 

Ultrasonic  Images 

Udantha  R.  Aheyratne,  Athina  P.  Petropulu* ,  John  M.  Reid  and  Thomas  Colas* 

Biomedical  Engineering  and  Science  Institute 
*Department  of  Electrical  and  Computer  Engineering, 

Drexel  University,  Philadelphia,  PA  19104 

We  address  the  problem  of  estimating  imaging-distortions  and  modeling  the  tissue  re¬ 
sponse  from  clinical  B-scan  ultrasound  images,  based  on  rf  A-lines  that  form  the  image.  We 
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model  rf  A-lines  as:  yi{n)  =  hi{ri)  *  fi{n)  +  Wi{n),i  =  1,2,  •  •  •,  where  yi(n)  is  the  observed 
rf  A-line,  fi{n)  is  the  tissue  response,  Wi{n)  is  the  observation  noise  and  /i,(n)  is  the  ker¬ 
nel  representing  imaging  distortions;  the  symbol  stands  for  the  convolution  operation. 
The  fol-  lowing  assumptions  are  made:  [Al]  hi{n)  is  deterministic,  possibly  non-  minimum 
phase,  [A2]  fi{n)  is  stationary,  zero-mean  non-Gaussian,  modeled  as  the  convolution  of  a 
white  component  ej(n)  and  a  deterministic  part  ti{n)  corresponding  to  the  statistical  color 
of  the  issue  response,  [A3]  Wi{n)  is  zero-mean,  Gaussian,  and  independent  of  fi{n).  A  col¬ 
ored  random  process  is  taken  as  the  model  for  the  tissue  response,  considering  the  fact  that 
under-  lying  scatterers  need  not  be  randomly  and  independently  located.  Applying  higher 
order  spectra  based  blind  deconvolution  on  pairs  of  closely  located  rf  A-lines,  we  identify 
the  color  of  the  tissue  response  as  well  as  the  distortion  kernels.  The  color  of  the  tissue 
response  is  independent  of  the  imaging  system  and  can  be  used  in  tissue  characterization. 
We  estimate  distortion  kernels  and  the  color  of  the  tissue  response,  from  a  B-mode  image 
of  a  human  liver,  containing  a  tumor.  The  statistical  color  of  the  tissue  response  estimated 
inside  the  tumor  has  features  that  are  distinctly  different  from  those  estimated  outside  the 
tumor.  This  may  be  attributed  to  the  change  in  the  tissue  structures  brought  about  by 
the  tumor.  Our  kernel  estimation  and  color  identification  process  is  based  on  data  from 
regions  as  small  as  12  mm^  (approx.)  in  the  tissue  space.  The  advantages  of  the  proposed 
approach  are: (a)  since  the  estimations  are  carried  out  in  higher  order  spectra  domain,  results 
are  immune  to  additive  Gaussian  noise  in  observations,  and  (b)  both  minimum  phase  and 
non- minimum  phase  kernels  can  be  estimated.  We  compare  our  method  with  existing  non- 
parametric  pulse  estimation  techniques  and  demonstrate  relative  strengths  and  weaknesses. 


IEEE  Workshop  on  Nonlinear  Signal  and  Image  Processing,  Neos  Marmaras,  Greece,  June 
1995. 


Higher  Order  Spectra  Based  Deconvolution  of  Ultrasound  Images 

Udantha  R.  Abeyratne,  Athina  P.  Petropulu'^,  John  M.  Reid 
Biomedical  Engineering  and  Science  Institute 
1  Electrical  and  Computer  Engineering  Department 
Drexel  University,  Philadelphia,  PA  19104 

Our  goal  is  to  model  and  identify  the  tissue  response  based  on  the  backscattered  rf 
signal  that  forms  the  ultrasound  image.  We  model  the  rf  ultrasonic  image,  in  both  axial 
and  lateral  directions,  as  the  convolution  between  an  1-D  hypothetical  tissue  response  and 
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the  ultrasonic  system  response  (distortion),  plus  additive  Gaussian  noise.  We  model  the 
tissue  response  as  a  non-Gaussian,  colored  random  process.  Closely  spaced  axial  (lateral) 
image  lines  contain  as  a  common  convolutional  term  the  axial  (lateral)  distortion,  whereas 
the  noncommon  terms  are  due  to  the  tissue  color.  Applying  blind  deconvolution,  we  identify 
the  color  of  the  tissue  response  as  well  as  the  corresponding  distortion  kernels.  The  color 
of  the  tissue  response  is  independent  of  the  imaging  system  and  can  be  used  as  a  tissue 
characterization  parameter.  Estimated  kernels  are  compared  with  experimentally  obtained 
kernels,  and  the  deconvolution  of  real  ultrasound  images  is  performed. 
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a- Weighted  Cumulant  Projection:  A  New  Tool  For  System  Identification 

Udantha  R.  Abeyratne  and  Athina  P.  Petropulu* 

Biomedical  Engineering  and  Science  Institute 
*Department  of  Electrical  and  Computer  Engineering, 

Drexel  University,  Philadelphia,  PA  19104 

Summary 

It  is  well  established  that  identification  of  a  nonminimum  phase  non-Gaussian  random 
process  can  be  achieved  based  on  its  higher-order  cumulants  (order  three  or  higher).  We 
propose  a  computationally  attractive  approach  for  system  identification,  whose  complex¬ 
ity  remains  constant  as  the  order  of  the  employed  cumulants  increases.  We  define  the 
a— weighted  n*'^-order  cumulant  projection  of  the  process  x{k)  to  be 

«)  =  •••  m  "^2,  •••,  ^  ^  .  complex,  |q;|  =  1,  (17) 

r2  Tn-l 

where  c®(-)  denotes  order  cumulant.  For  the  process 

x{k)  =  h{k)  *  e{k),  (18) 

where  e{k)  is  white  non-Gaussian  zero-mean,  and  h{k)  is  deterministic  generally  nonmini¬ 
mum  phase,  we  prove  that 

P^{z-,  a)  =  (19) 

where  P^{z\a)  is  the  Z-transform  of  p®(r;a),  and  c  is  a  constant.  ^From  (19),  it  can  be 
shown  that  by  controlling  a,  H{z)  can  be  identified  for  n  >  2.  If  p^(A;;  a)  and  h{k)  denote 
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the  complex  cepstra  of  p®(r;  a)  and  h{k),  respectively,  we  show  that 


(20) 

(21) 

which  leads  to  a  scaled  and  shifted  version  of  h{k)  via  inverse  cepstrum  operations.  We 
show  that  we  can  always  find  a  complex  number  a  with  unit  magnitude  to  guarantee  that 
the  denominator  in  (20)  and  (21)  is  different  from  zero,  and  discuss  the  rationale  behind 
choosing  a  particular  value  for  a. 

Based  on  equation  (19)  a  formula  for  the  estimation  of  the  Fourier  phase  of  h{k)  from 
the  phase  of  the  a- weighted  cumulant  projection  is  also  derived. 
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Higher-Order  Statistics  Based  Approaches  for  Estimating  Scatterer 
Periodicity  and  Correlation  Structure  of  Tissue 

Udantha  R.  Abeyratne,  Athina  P.  Petropulu* ,  John  M.  Reid,  Ethan  J.  Halpern'^  and 

Flemming  Forsherg'^ 

Biomedical  Engineering  and  Science  Institute 
*  Electrical  and  Computer  Engineering  Department  Drexel  University,  Philadelphia,  PA 

19104 

+Dept.  of  Radiology,  Thomas  Jefferson  University  Hospital,  Philadelphia.  PA  19107. 

We  model  tissue  as  a  collection  of  point  scatterers  embedded  in  a  uniform  media,  and 
show  that  the  higher-  order  statistics  of  the  scatterer  spacing  distribution  can  be  estimated 
from  digitized  RF  scan  line  segments  and  be  used  in  obtaining  tissue  signatures.  We  model 
the  ultrasonic  RF  echo,  y(t),  by  y(t)=h(t)*f(t)-(-w(t),  where  h(t)  is  the  pulse-echo  wavelet 
of  the  imaging  system,  w(t)  is  the  zero-mean  observation  noise,  and  is  the  linear  con¬ 
volution  operator.  The  quantity  f(t),  which  we  will  call  the  tissue  response,  represents 
the  underlying  tissue  structure,  and  is  modeled  as:  f(t)  =  rd(t)-frm(t)-l-rr(t),  where,  rd(t) 
represents  unresolvable  diffuse  scatterers  leading  to  fully-developed  speckle  with  Gaus¬ 
sian  statistics  [1],  [2],  rr(t)  represents  resolvable,  coherent  scatterers  with  long-range  order 
(pseudo-periodicity). The  quantity  rm(t)  represents  unresolvable  periodicity  from  structural 
scatterers,  and,  correlated  non-periodic  components  of  diffused  and  structural  scatterers 
(short-  range  order).  The  terms  rm(t),  rr(t)  (and  their  contributions  to  the  RF-echo  y(t)) 
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are  non-Gaussian,  due  to  contributions  from  periodic  structures  and/or  scatterers  of  short 
range  order.  Based  on  our  model  for  tissue  micro  structure,  we  develop  methods  for  the 
estimation  of  resolvable  periodicity  via  higher-order  spectra  of  the  RF  observation  y(t). 
Also,  we  propose  a  tissue  signature,  the  color  of  the  tissue  response  ,  that  captures  the 
correlation  among  non-periodic  scatterers.  The  tools  employed,  i.e.,  higher-order  statis¬ 
tics,  were  chosen  as  the  most  appropriate  ones  because  they  suppress  Gaussian  processes, 
such  as  the  one  arising  from  the  diffused  scatterers  and  observation  noise.  Higher-order 
statistics,  unlike  second-order  statistics,  also  preserve  the  Fourier-phase  of  the  color  of  the 
tissue  response.  Working  on  simulated  and  clinical  data,  we  show  that  the  proposed  peri¬ 
odicity  estimation  technique  is  superior  to  the  widely  used  power  spectrum  and  cepstrum 
techniques  in  terms  of  the  accuracy  of  estimations.  We  also  show  that  even  when  there  is 
no  significant  periodicity  in  data,  we  are  still  able  to  characterize  tissues  using  signatures 
based  on  the  higher-order  cumulant  structure  of  the  scatterer  spacing  distribution. 

[1]  K.D.  Donohue,  J.M.  Dressier,  T.  Varghese  and  N.M.  Bilgutay,  “Spectral  Correlation 
in  Ultrasonic  Pulse-echo  Signal  Processing”,  IEEE  Tr.  on  Ultrasonics,  Ferroelec.,  and 
Frequency  Control,  vol.  40,  no.  4,  pp. 330-337,  1993. 
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On  Modeling  The  Tissue  Response  From 
Ultrasonic  B-scan  Images 

Udantha  R.  Abeyratne,  Athina  P.  Petropulu,  John  M.  Reid 


Abstract —  We  model  tissue  as  a  collection  of  point  scat- 
terers  embedded  in  a  uniform  media,  and  show  that  the 
higher-order  statistics  {HOS)  of  the  scatterer  spacing  dis¬ 
tribution  can  be  estimated  from  digitized  RF  scan  line  seg¬ 
ments  and  be  used  in  obtaining  tissue  signatures.  We  as¬ 
sume  that  RF  echoes  are  non-Gaussian,  on  the  grounds  of 
empirical/theoretical  justifications  presented  in  works  such 
as  [25],  [6],  [3],  [24]  and  [27].  Based  on  our  model  for  tis¬ 
sue  microstructure,  we  develop  schemes  for  the  estimation 
of  resolvable  periodicity  as  well  as  correlations  among  non¬ 
periodic  scatterers.  Using  higher-order  statistics  of  the  scat¬ 
tered  signal,  we  define  as  tissue  ‘‘color”  a  quantity  that  de¬ 
scribes  the  scatterer  spatial  correlations,  show  how  to  evalu¬ 
ate  it  from  the  higher-order  correlations  of  the  digitized  RF 
scan  line  segments,  and  investigate  its  potential  as  a  tissue 
signature.  The  tools  employed,  i.e.,  higher-order  statistics, 
were  chosen  as  the  most  appropriate  ones  because  they  sup¬ 
press  Gaussian  processes,  such  as  the  one  arising  from  the 
diffused  scatterers.  Higher-order  statistics,  unlike  second- 
order  statistics,  also  preserve  the  Fourier-phase  of  the  signa¬ 
ture,  the  color  of  the  tissue  response.  Working  on  simulated 
and  clinical  data,  we  show  that  the  proposed  periodicity 
estimation  technique  is  superior  to  the  widely  used  power 
spectrum  and  cepstrum  techniques  in  terms  of  the  accuracy 
of  estimations.  We  also  show  that  even  when  there  is  no  sig¬ 
nificant  periodicity  in  data,  we  are  still  able  to  characterize 
tissues  using  signatures  based  on  the  higher-order  cumulant 
structure  of  the  scatterer  spacing  distribution. 

Keywords — tissue  modeling,  tissue  characterization,  med¬ 
ical  ultrasound,  higher-order  statistics,  periodicity  estima¬ 
tion 


1.  Introduction 

Ultrasound  is  one  of  the  widely  used  medical  imag¬ 
ing  techniques,  due  to  its  versatility,  relative  safety, 
low  cost  and  the  availability  of  portable  units.  It  is  com¬ 
monly  used  by  radiologists  as  an  adjunct  modality  in  dif¬ 
ferentiating  between  cystic  and  solid  masses,  a  crucial  dis¬ 
tinction  in  diagnosing  cancer.  Despite  decades  of  research 
leading  to  many  technological  advances,  still  the  clinical  di¬ 
agnoses  are  primarily  based  on  qualitative  evaluation  of  im¬ 
ages  by  expert  observers.  There  have  been  many  attempts 
at  developing  objective  tissue  characterization  criteria  on 
the  premise  that  there  is  much  more  observer-independent 
information  available  from  ultrasound  than  what  is  cur- 
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rently  being  used.  These  are  rooted  on  the  fundamental 
notion  that  the  biological  tissues  are  composed  of  charac¬ 
teristic  structures  whose  ultrasonic  properties  often  change 
due  to  diseases.  The  goal  of  tissue  characterization  is  the 
extraction  of  signatures  that  assume  distinct  values  in  the 
presense  of  normal  and  diseased  states  of  tissues,  such  that 
it  is  possible  to  differentiate  between  them. 

Although  the  exact  identities  of  the  physical  structures 
responsible  for  ultrasound  backscattering  are  generally  un¬ 
known  [22],  [5],  they  can  usually  be  characterized  either 
as  macroscopic  or  microscopic,  compared  with  the  wave¬ 
length  of  the  interrogating  ultrasonic  beam.  The  tissue 
is  often  modeled  as  a  collection  of  point  scatterers,  em¬ 
bedded  in  a  uniform  non-scattering  medium.  Considering 
the  biological  variabilities  associated  with  tissues,  the  spa¬ 
tial  distribution  and  the  scattering  strengths  (scattering 
cross-sections)  associated  with  these  scatterers  are  usually 
described  in  statistical  terms.  The  statistics  of  the  inter- 
scatterer  spacing  distribution  are  commonly  used  as  tissue 
signatures. 

In  organs  such  as  the  liver  the  overall  arrangement  of 
the  tissues  consists  of  an  organized  repetition  of  a  basic 
structural  unit.  This  repeated  structure  is  regarded  to  be 
providing  resolvable,  repeated  scattering  centers  for  the 
propagating  ultrasonic  pulse  [10],  [12],  [13],  [7].  It  has 
been  shown  that  a  periodicity  can  be  observed  from  the 
liver  pulse-echo  data,  and  the  estimated  periodicity  can 
be  linked  to  the  mean-scatterer-spacing  of  the  underlying 
tissue  [10].  It  was  also  shown  that  the  mean  scatterer  spac¬ 
ing  may  be  used  as  a  feature  for  tissue  characterization  of 
diffuse  liver  diseases,  such  as  cirrhosis  and  chronic  active 
hepatitis  [10], [13]. 

However,  the  success  of  some  of  these  methods  [10]  was 
reported  to  be  compromised  by  the  presence  of  coherently 
reflecting  structures  such  as  blood  vessels.  In  [26],  the  pe¬ 
riodicity  detection  was  carried  out  in  the  power  cepstrum 
domain.  In  [13],  the  liver  tissue  micro-structure  has  been 
described  by  a  three  component  model:  the  first  compo¬ 
nent  (diffuse)  represents  randomly  positioned  scatterers 
of  sufficient  concentrations  to  generate  circular- Gaussian 
statistics;  the  second  component  (regular)  represents  non- 
randomly  distributed  coherent  scatterers  with  long-range 
order  (pseudo-periodicity);  the  third  component  (short- 
range  order)  is  due  to  organ  boundaries  and  blood  vessels 
etc. 

All  of  the  above  techniques  are  based  on  the  auto- 
correlation/power  spectrum  of  the  observations,  thereby 
neglect  information  contained  in  the  Fourier-phase  of  the 
RF  echo.  In  [28],  [9]  the  RF  echo  was  modeled  as  a  cyclo¬ 
stationary  signal,  thus  preserving  phase  information.  The 


IEEE  TRANSACTIONS  ON  MEDICAL  IMAGING,  AUGUST  1996 

embedded  periodicity  was  detected  in  the  frequency  do¬ 
main.  That  approach  seemed  to  perform  better  than  con¬ 
ventional  power-spectral  density/cepstra  based  techniques 
in  detecting  periodicity.  Both  in  the  cyclo-stationarity  ap¬ 
proach  and  the  power  cepstrum  one  it  was  assumed  that 
the  diffused  scatterer  spacings  are  un correlated. 

The  scatterer  spatial  correlation  has  been  used  in 
the  past  to  model  ultrasonic  properties  of  tissue 
[11], [15], [16], [8], [4]  and  was  shown  to  lead  to  promising  tis¬ 
sue  signatures.  In  [16],  it  was  shown  that  using  power 
spectra,  deterministic,  membrane-like  structures  as  well 
as  structures  consisting  of  random,  diffuse  structures  can 
be  characterized.  In  [17],  effective  scatterer  sizes,  concen¬ 
trations  and  acoustic  impedances  were  investigated  using 
power  spectra,  as  potential  tissue  signatures.  The  effective 
scatterer  size  was  reported  to  be  the  most  important  tissue 
feature  sensed  with  the  method  of  [17]. 

In  this  paper,  we  use  a  point  scatterer  model  similar  to 
the  ones  in  [13]  and  [9]  to  describe  tissue  microstructure, 
with  the  exception  that  our  model  takes  the  correlations 
existing  among  both  resolvable  and  non-resolvable  scatter- 
ers  into  account.  This  enables  us  to  also  consider  cases 
of  coherent  long-range  scatterer  distributions  which  have 
high  variances  associated  with  the  inter-scatterer  spacing, 
i.e.  almost  no  periodicity,  and,  non-periodic  echos  resulting 
from  correlated,  non-periodic  components  of  both  diffused 
and  structural  scatterers.  We  assume  that  RF  echoes  are 
non-Gaussian,  on  the  grounds  of  empirical/theoretical  jus¬ 
tifications  presented  in  works  such  as  [25],  [6],  [3],  [24]  and 
[27].  Based  on  our  model  for  tissue  microstructure,  we  de¬ 
velop  schemes  for  the  estimation  of  resolvable  periodicity  as 
well  as  correlations  among  non-periodic  scatterers.  Using 
higher-order  statistics  of  the  scattered  signal,  we  define  as 
tissue  “color”  a  quantity  that  describes  the  scatterer  spatial 
correlations,  show  how  to  evaluate  it  from  the  higher-order 
correlations  of  the  digitized  RF  scan  line  segments,  and 
investigate  its  potential  as  a  tissue  signature.  The  tools 
employed,  i.e.,  higher-order  statistics,  were  chosen  as  the 
most  appropriate  ones  because  they  suppress  Gaussian  pro¬ 
cesses,  such  as  the  one  arising  from  the  diffused  scatterers. 
Higher-order  statistics,  unlike  second-order  statistics,  also 
preserve  the  Fourier-phase  of  the  signature,  the  color  of  the 
tissue  response. 

Via  computer  simulations,  we  show  that  the  proposed 
periodicity  estimation  technique  is  superior  to  the  widely 
used  power  spectrum  and  cepstrum  techniques  in  terms 
of  the  accuracy  of  estimations.  We  also  preliminary  evi¬ 
dence  that  even  when  there  is  no  significant  periodicity  in 
data,  we  are  still  able  to  characterize  tissues  using  signa¬ 
tures  based  on  the  higher-order  cumulant  structure  of  the 
scatterer  spacing  distribution. 

II.  Theory 

In  the  following,  we  use  lower  case  characters  to  symbol¬ 
ize  time  domain  quantities,  and  upper  case  ones  to  denote 
frequency  domain  quantities  unless  otherwise  specifically 
specified. 

Assuming  a  narrow  ultrasound  beam,  linear  propagation 


2 

and  weak  scattering  in  the  tissue  medium,  we  model  the 
ultrasonic  RF  echo,  y(f),  by, 

y{t)  =  h{t)*  f{t)  +  w{t),  (1) 

where  h(t)  is  the  pulse-echo  wavelet  of  the  imaging  sys¬ 
tem,  w{t)  is  the  zero-mean  observation  noise  and  is  the 
linear  convolution  operator.  The  quantity  f(t)  represents 
the  underlying  tissue  structure  and  will  be  called  the  tissue 
response. 

A.  The  tissue  structure: 

Following  [5], [13]  and  [9]  we  model  the  structures  within 
tissue  that  are  responsible  for  backscattered  ultrasonic  field 
by  point-scatterers  organized  at  different  levels.  Existing 
work  such  as  [13]  and  [9]  attempt  to  characterize  liver  tis¬ 
sue  based  on  a  two  component  point  scatter  model:  (a) 
the  diffused  component^  which  represents  [13]  randomly  po¬ 
sitioned  scatterers  of  sufficient  concentrations  to  produce 
echo  signals  with  circular  Gaussian  statistics.  The  po¬ 
sitions  of  individual  scatterers  are  assumed  uncorrelated 
[13];  (b)  the  coherent  component^  which  represents  non- 
randomly  distributed  scatterers  with  long-range  order  [13], 

In  this  paper,  we  describe  liver  tissue  by  a  three- 
component  model,  which  in  addition  to  (a)  and  (b)  above, 
also  takes  into  account  the  correlations  existing  among  (i) 
diffused  scatterer  locations,  (ii)  resolvable  and  unresolvable 
pseudo-periodic  scatterers  leading  to  long  range  order,  and 
(iii)  collection  of  scatterers  leading  to  short  range  order, 
but  not  strong  enough  to  violate  the  conditions  of  weak- 
scattering  or  the  stationarity  of  the  RF  echo  in  a  significant 
manner. 

We  model  the  tissue  response  f{t)  by, 

fit)  =  rdit)  +  rmit)  +  rr{t),  (2) 

where, 

•  rd{t)  =  (ip6{t  —  i/p)  represents  unresolvable  dif¬ 

fuse  scatterers  leading  to  fully-developed  speckle  with 
Gaussian  statistics, 

•  rm(t)  =  ^p)  represents  the  combined  ef¬ 

fects  of  unresolvable  periodicity  from  structural  scat¬ 
terers,  and,  correlated  non-periodic  components  of  dif¬ 
fused  and  structural  scatterers. 

•  the  term  rr{t)  =  Ylp~i  ~  ^p)  represents  the  peri¬ 
odic,  resolvable  scatterers. 

Note  that  the  terms  rm{i)  and  rr{t)  are  non-Gaussian,  due 
to  contributions  from  periodic  structures  and/or  scatterers 
of  short  range  order. 

B.  The  RF-echo: 

From  (1)  and  (2),  y{t)  can  be  expressed  as: 

y{i)  ■=  yd{t)  +  ym{t)  +  yr{i)  +  (3) 

where  yd{t)  =  rd{t)^h{t),  ym{t)  =  rm{t)^h{t)  and  yrit)  = 
rr{i)  *  h{t)  are  respectively  called  the  diffused,  mixed  and 
resolvable-periodic  components  of  the  RF  echo  y{t).  The 
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observation  noise  w{t)  is  assumed  zero-mean,  Gaussian  and 
uncorrelated  with  yi{t),  i  = 

Let  us  assume  an  infinite  extent  tissue  media,  and 
Nd^Nm^Nr  oo.  The  RF  echo  y{t)  is  assumed  to  be 
compensated  for  attenuation.  In  section  2.3,  we  discuss 
the  effects  of  attenuation  and  finite-extent  data  on  our  es¬ 
timations. 

Assuming  that  the  terms  yi{t)j  i  =  d^m^r  are  zero- 
mean,  mutually  uncorrelated,  the  third-order  cumulants 
C3(Ti,r2)  [19]  the  process  y{t)  can  be  written  as: 

cI{ti,T2)  -  E  {y{t)y{t  +  Ti)y{t  +  T2)}  , 

=  4‘'(n,T2)  +  c|’”(ri,r2) 

+4"(n,'r2),  (4) 

where  E{’}  is  the  statistical  expectation  operator  over  the 
ensemble.  The  third-order  cumulants  Cg  (ri ,  T2)  of  the  noise 
process  w{t)  vanishes  under  the  hypothesis  that  w{t)  is 
Gaussian. 

B.l  The  diffused  component  (DC): 

We  model  the  diffused  component  yd{t)  of  the  RF  echo 
y{t)  as  a  Gaussian  process  (See  also  [13]),  As  such  the 
third-order  cumulants  sequence  equals: 

C3‘'(ri,r2)  =  0,  Vri,r2.  (5) 

B.2  The  mixed,  non-Gaussian  component  (MC): 

The  mixed  component  i/m(0  of  the  RF  echo  y{t)  is  given 
by: 


Let  us  assume: 

(Al)  Xi  is  uncorrelated  to  bj,  i,j  =  1,  2,  •  •  • ,  Am, 

(A2)  Xi^  i  =  1, 2,  •  •  form  a  stationary  process,  and, 
(A3)  2/m(f)  is  a  zero-mean  process. 

Under  (A3)  the  third  order  cumulants  of  ym{i)  equal, 

cl^Cn,  T2)  =  E{ymit)ym{t  +  ri)ym{t  +  T2)], 

=  M3"‘(ri,r2)*c|(ri,r2),  (7) 

with, 

Nrr.  Nm 

M;-iruT2)  = 

p=:l  q  =  l  s  =  l 

xE{6(t  -  Xp)6{t  +  ri  -  Xq)6{t  +  r2  -  A^)}  ,  (8) 

where  (Al)  has  been  used  in  obtaining  (8). 

Let  the  joint  probability  density  function  (arbitrary) 
for  the  scatterer  location  triplets  Xp^Xg^Xs  be  given  by 
/m(Ap,  A^,  A5).  Then  from  (8)  we  obtain: 

Nm  Nm 

M^^(rur2)  = 

p=l  q=l  s=l 


00  00  00 

^/// 


—  Ap)^(i  Ti  —  Xq)6{t  +  T2  —  A5 


-00  —00  —00 

X/m(Ap,  Agi  X  s)d\pdXqd\$ 

Nm  Nm  Nm 

P=1  q=zl  S  =  1 

+  T2). 


(9) 


Under  the  assumption  (A2),  the  probability  density  func¬ 
tion  /m(fjf  +  +  T2)  does  not  depend  on  the  particu¬ 
lar  t  chosen.  Thus,  /m(f,^  +  +  ^2)  can  be  written  as 

/m('3^i,  7’2)-  Therefore,  from  (7)  and  (9)  we  obtain: 

C3’"(ri,r2)  =  Kmfm{Tl,r2)  *  C3{ti,T2),  (10) 

where  Km  is  given  by, 

Nm  Nm  Nm 

^-  =  EEE^{WJ-  (11) 

p=l  q=l  r=:s 


B.3  Resolvable  periodic  component  (RPC): 

In  the  context  of  ultrasonic  tissue  characterization,  the 
RPC  component  has  been  shown  to  be  of  great  significance 
[10],  [9].  To  be  classified  as  the  resolvable  periodic  com¬ 
ponent,  the  scatterer  separation  should  be  regular  enough, 
and  the  repeat  distances  should  be  large  enough  (compared 
to  the  length  of  the  pulse-echo  impulse  response)  to  he  re¬ 
solvable.  The  Gamma  distribution  has  proven  to  be  a  par¬ 
ticularly  useful  tool  in  describing  inter-scatterer  space  dis¬ 
tribution  [14],  [9],  because  of  its  versatility  in  producing 
scatterer  locations  ranging  from  almost  periodic  to  clearly 
non-periodic. 

The  general  mathematical  description  of  the  RPC  fol¬ 
lows  directly  from  the  expressions  derived  for  the  case  of 
MC  component.  The  resolvable  periodic  component,  yr{t)^ 
can  be  written  as: 


Nr 


—  Op) 


*  h{t). 


(12) 


Again,  we  assume: 

(Bl)  Oi  is  uncorrelated  to  vj,  i,j  =  1, 2,  •  •  • ,  , 

(B2)  Oi,  z  =  1, 2,  •  •  froms  a  stationary  process,  and, 
(B3)  yr(f)  is  a  zero-mean  process. 

From  (12)  we  can  obtain  as  in  the  case  of  MC  compo¬ 
nent, 


C3’'(^1.^2)  =  Krfr{Tl,T2)*C2{ri,T2),  (13) 

Nr  Nr  Nr 

Kr  -  EEE-^(^P^?^‘'L  (1^) 

p=l  qzzl  5  =  1 


where  the  joint  probability  density  function  /r(ri,  r2)  now 
describes  a  (pseudo)  periodic  phenomenon,  i.e.  /r('ri,  ^2)  is 
a  two  dimensional  bed  of  spikes,  whose  separation  is  equal 
to  the  mean-scatterer-spacing. 
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In  the  following  we  derive  an  expression  for  fr{Ti,r2)^ 
under  the  constraint  that  the  resolvable  periodic  scatter- 
ers  are  separated  by  a  constant  time  interval,  T,  i.e.  the 
process  is  strictly  periodic  and  the  resolvable  periodicity 
(unknown)  is  T.  Then  the  periodic  component  rr{t)  can 
be  written  as: 

Nr 

=  (15) 

pz=l 

where  Vp{t)  is  the  scaiterer  strength  function^  and  is  a 
bounded  function  defined  over  the  set  of  real  numbers, 

^(*“00,0o)  • 

The  Fourier  transform  Rr{Lo)  of  rr{t)  can  be  obtained  as, 

=  ^  f]  V{L,-^k),  (16) 

k=—oo 

where  Vp{Lo)  is  the  Fourier  transform  of  Vp{t).  Let  the 
region  of  support  of  the  function  Vp{t)  be  {t  :  \t\  <  B},  As 
5  — )■  oo,  V{oj)  — ^  h{uj).  As  B  ^  oo  the  bispectrum  [19] 
C^'"{u}i,U2)  of  rr{i)  is  given  by: 

=  Rr{uJi)Rr{i02)Rr{—^l—<^2)y 

ki  k2  ks 

27r 

x6{-W2  -  Wi  -  — fcs),  (17) 

=  ( ^  ~  Y^2)-  (18) 

ki  k2 

For  the  case  V{io)  — ^  ^((^)>  from  (12)  and  (18)  we  can 
obtain  the  bispectrum  Cl''{oJiyU)2)  of  yr{t)  as: 

Ct{wi,W2)  = 

ki  k2 

2'7r 

-  ^^2)}-  (19) 

The  third-order  cumulant  sequence  C3'‘(ri,T'2)  of  the 
RPC  can  be  obtained  from  (19)  through  an  inverse  Fourier 
transform,  i.e., 

<?z''{ti,T2)  =  4(n,T2)  -  k{r)8{T2  -  k2T). 

ki  k2 

(20) 

Comparing  (14)  and  (20)  we  get: 

fr{Tl,T2)  =  Y,Y.^{ri  -  klT)6{T2  -  k2T),  (21) 

A;i  k2 

Based  on  (5),  (10),  and  (13),  eq.  (4)  becomes: 

C3  (ti  ,  T2  )  =  {Km  fm  (ti  ,  r2  )  +  AT  /r  (ti  ,  T2  )  }  *  C J  (ti  ,  T2  ) , 

=  {Kmfm{ri,T2)  +  YjY,  ■  ^i7’)«(t2  -  hT)} 

ki  k2 

*C3(ri,r2),  (22) 

when  the  inter-scatterer  spacing  of  the  RPC  component 
is  a  constant  T. 


(7.  The  proposed  method  for  the  estimation  of  periodicity 

From  (18)  and  (22)  we  observe  that  the  periodicity  man¬ 
ifests  itself  in  both  bispectrum  and  the  third  order  cumu¬ 
lant  domains.  In  the  cumulant  domain,  peaks  are  separated 
by  the  period  T,  while  in  the  bispectrum  domain  the  sep¬ 
aration  is  ^  rad/S.  Let  us  form  the  weighted  S’^^-order 
cumulant  projection  (see  also  [2],[1])  of  (22)  as  follows: 

pK^.w)  =  ^c|(rl,r2)e■'“’■^  0  <  w  <  27r,  (23) 

7-2 

=  (24) 

From  (23),  pKri,^;)  can  be  viewed  as  the  one-dimensional 
Fourier  transform  of  c^ri,  T2)  with  respect  to  the  variable 
T2.  According  to  this  interpretation,  p^(ri,a;)  equals: 

Ps{ri,uj)  =  Fr2{cl"'{T’i>r2)] , 

=  Ayn-Fra  {/m(Tl ,  T2)  *  C3  (ti,  r2)}  ,  (25) 

and  also 

vl{ruuj)  =  Fr2  {(c^"(n,T2))}  , 

=  ^fc2),  (26) 

fci  k2 

where  FrA-}  is  the  one-dimensional  Fourier  transform 
with  respect  to  the  variable  7-2,  and  C3”^(ri,r2)  and 
C3''(ri,r2)  are  given  by  (10)  and  (20). 

Based  on  (26)  the  periodicity  can  be  estimated  from  the 
p|(ri,a;)  as  the  separation  T  between  peaks  along  the  ti~ 
axis.  Along  the  w-axis,  the  peaks  are  separated  by  the 
amount  27: fT, 

D.  The  estimation  of  the  color  of  the  tissue  response 

When  the  tissue  under  study  shows  no  periodic  varia¬ 
tions  in  ultrasonic  properties,  the  periodicity  will  loose  its 
significance  as  a  tissue  signature.  In  the  following  we  define 
a  new  quantity,  the  color  of  the  tissue  response  to  supple¬ 
ment  periodicity  as  a  tissue  signature. 

The  contribution  of  the  MC  to  the  third-order  cumu- 
lants  of  the  RF  echo  is  described  by  (10).  The  term 
Kmfm{T'i^T’2)  Can  be  viewed  as  the  third-order  moment  se¬ 
quence  of  the  signal: 

s{t)  =  e(t)  *g{t),  (27) 

where  e{t)  is  an  i.i,d.^  non-Gaussian,  zero-mean  noise  pro¬ 
cess,  whose  skewness  is  Km’,  and,  g{t)  is  the  signal  whose 
third-order  moment  sequence  is  given  by  /m(Ti,r2). 

We  will  refer  to  ^(if)  as  the  “color  of  the  tissue  response” . 
The  time  domain  equivalent  of  (10)  becomes: 

z{t)  =  e{t)  *  u{t)j  (28) 

where  u{t)  =  g{t)  ^h{t)  will  be  referred  to  as  the  combined 
kernel 

Assuming,  that  the  support  of  the  third-order  cumulant 
sequence  C3(ri,r2)  of  z(t)  is  contained  in  a  square  defined 
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by  Ini,  |r2|  <  T,  where  T  is  the  resolvable  periodicity,  we 
can  recover  the  term  c^rijn))  from  cKn^n)  as: 

cI{ti,T2)  =  c|(ri,r2),  Ini,  Ini  <  Tq,  (29) 

where  Tq  ^  T  and  T  is  the  known  or  estimated  period  of 
the  RF  data.  Note  that  a  precise  knowledge  of  T  is  not 
required. 

Applying  methods  of  system  reconstruction  from  higher- 
order  spectra,  [19],  [2],  [1],  we  can  retrieve  u{t)  from 
cKnjn)  within  a  scalar  and  a  time  delay.  In  the  work 
of  this  paper,  we  used  the  a-WCP  based  method  of  system 
reconstruction  proposed  in  [2],  exploiting  its  low  computa¬ 
tional  load. 

The  color  of  the  tissue  response,  g{t)j  compactly  sum¬ 
marizes  both  phase  and  magnitude  information  available 
in  the  tissue  response  f{t). 

The  quantities  Ui{i)  and  Uj{t)  obtained  from  two  differ¬ 
ent  A-lines  and  “j”  can  be  written  as: 

Ui{t)  =  h{t)*gi{t),  (30) 

Uj{t)  =  h{t)*gj{t),  (31) 

where  h{t),  the  pulse-echo  impulse  response  of  the  imaging 
system,  can  be  assumed  to  be  invariant  between  the  two 
locations  ‘‘i”  and  gj{t)  respectively  represent 

the  color  of  the  tissue  responses  at  “z”  and  .  The  quan¬ 
tities  gk{i),  k  =  1,2,  summarizes  information  about 
the  tissue,  and  is  also  free  from  the  effects  of  making 
it  a  good  candidate  for  a  tissue  signature. 

Using  (30)  and  (31)  we  form  a  tissue  signature  Sij{t)^  in 
the  time  domain,  as  follows: 

=  h{t)*gi{t)*{h{t)*gj{t)}~^, 

=  9i{t)*g]'^it),  (32) 

where  {•}“^  denotes  the  convolutional  inverse. 

If  the  estimations  are  from  two  regions  with  similar  ul¬ 
trasonic  properties,  for  instance  from  two  healthy  regions 
of  a  liver,  Sij{t)  will  approach  a  Dirac  delta-function.  On 
the  otherhand,  if  the  two  estimates  are  from  two  regions  in¬ 
side  and  outside  a  tumor  which  has  brought  in  a  change  in 
the  inter-scatter  spacing  distribution,  then  Sij{t)  and  will 
depart  from  Dirac  delta.  Thus,  the  closeness  of  Sij  to  a 
Dirac  delta  performs  as  a  relative  tissue  signature.  We  will 
estimate  Sij  from  clinical  images  and  validate  its  usefulness 
as  a  signature. 

The  quantities  gi{t)  and  gi{t)  can  individually  be  ob¬ 
tained  from  (30)  and  (31)  using  a  blind  deconvolution  tech¬ 
nique  [21],  i.e.  a  technique  through  which  the  input  to  a 
system  can  be  obtained  when  only  its  output  is  known. 
Preliminary  results  of  a  blind  deconvolution  approach  to 
separate  gi{t)  and  gi{t)  have  been  reported  in  [4]. 

E.  The  effects  of  attenuation  and  finite  data  lengths 

The  results  in  Section  2  were  derived  under  the  assump¬ 
tion  of  an  infinite  tissue  media.  However,  in  practice  we 
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have  access  to  only  finite-length  data  records.  With  such 
records,  one  can  obtain,  parallel  to  the  situation  in  auto¬ 
correlation  estimation  from  finite  data,  only  a  statistical 
estimate  of  the  true  higher-order  cumulants  of  RF  data. 
Expressions  for  the  statistical  properties  of  HOS  estima¬ 
tors  can  be  found  in  [19].  Similar  to  the  power  spectrum 
domain,  in  the  bispectrum  domain  (and  in  pKti,^)),  the 
finite  data  records  limit  the  smallest  period  that  can  be 
detected. 

Attenuation  represents  all  forms  of  energy  loss  experi¬ 
enced  by  a  propagating  pulse  and  is  unavoidable  in  bio¬ 
logical  media.  The  effects  of  attenuation  on  the  pulse  at 
the  depth  x  —  where  v  is  the  ultrasonic  propagation 

velocity  and  to  is  time,  in  tissue  space  can  be  represented 
as  [18]: 

H{to,oj)=A{to,o^)H{io),  (33) 

where  H{to^uj)  and  A(^Oj^)  are  respectively  the  Fourier 
transforms  of  space-dependent  pulse  echo  response  and  at¬ 
tenuation;  H{lo)  is  the  Fourier  transform  of  the  pulse  trans¬ 
mitted. 

Eq.(33)  can  be  written  in  the  time  domain  as: 

h{to,t)  —  a{to,t)  *  h{t)y  (34) 

where  ^(^o,^),  a{to,t)  and  h{u))  are  the  time  domain  equiv¬ 
alents  of  A{toyt^)  and  H{u))  respectively. 

For  short  data  records,  we  can  ignore  the  to- dependency 
of  (34)  and  consider  as  a  time-invariant  quantity. 

Then  the  effect  of  attenuation  on  our  period  estimation 
techniques  follows  directly  from  (22),  and  (26)  with  the 
quantity  C3(ri,r2)  replaced  by  the  third-order  cumulant 
sequence  C3(ri,r2)  of  h(to,t)^ 

In  estimating  the  color  of  the  tissue  response,  the  com¬ 
bined  kernels  Ui{i)  and  Uj{t)  (see  (30)  and  (31))  now  takes 
the  form, 

Ui{t)  =  ai{tQ,t)^gi{t)  *  h{t),  (35) 

Uj{t)  =  aj{to)t)  *  gj{t)  *  h{t)^  (36) 

resulting  in  a  tissue  signature  Sij  given  by, 

Sij{t)  =  Ui(t)*{uj(t)}~^, 

=  (37) 

From  (37)  we  note  that  the  differences  in  attenuation 
characteristics  at  locations  “i”  and  “j”  are  also  captured 
in  Sij. 

III.  Data  collection 
A.  Simulated  data 

We  simulated  digitized  ultrasound  RF  scan  lines  on  a 
computer,  using  a  model  similar  to  the  ones  used  in  [14], 
[9]  and  [7].  Scatterers  were  modeled  as  a  collection  of 
points,  distributed  according  to  some  probability  density 
function.  A  uniform  probability  density  function  described 
both  inter-scatterer  spacings  and  the  scattering  strengths 
related  to  the  DC  component.  The  RPC  scatterer  spac¬ 
ings  were  described  by  a  F-probability  density  functions. 
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because  of  the  relative  ease  with  which  one  can  obtain  inter- 
scatterer  spacings  ranging  form  almost  periodic  to  uniform. 
Scatterer  strengths  were  also  modeled  with  a  F-density 
function  with  a  different  set  of  parameters.  The  quantity 
where  ^  is  the  mean  and  a  is  the  standard  deviation,  of  the 
scatter-strength  distribution  was  set  to  0.35,  so  that  we  get 
a  wide  variation  in  strengths.  The  program  also  had  the  ca¬ 
pacity  to  simulate  correlated  scatteres  contributing  to  the 
MC  component,  via  another  F-distribution.  The  correla¬ 
tions  in  the  scatterer  spacings  were  obtained  by  filtering 
them  through  a  linear,  time  invariant  system  of  arbitrary 
frequency  response. 

Observation  noise,  which  is  unavoidable  in  any  physical 
measurement,  was  modeled  by  a  zero-mean  Gaussian  pro¬ 
cess  as  is  commonly  done.  Observation  noise  was  allowed 
to  be  a  colored  process. 

The  center  frequency  of  the  transducer  was  A,2MHz] 
the  pulse  envelope  followed  a  Gaussian  shape  and  its  ZdB 
bandwidth  was  2,bMHz.  The  sampling  frequency  of  the 
data  acquisition  was  selected  to  be  20MHz. 

B.  Clinical  data 

To  demonstrate  the  performance  of  our  method  under 
clinical  conditions,  we  used  liver  scans  of  healthy  volun¬ 
teers,  (volunteer  nos.  I  and  II)  and  a  patient  (hypoe- 
choic  tumor  of  the  liver)  imaged  by  ultrasonologists  at 
the  Thomas  Jefferson  University  Hospital,  Philadelphia,  on 
a  model  UltraMark-9  system  manufactured  by  Advanced 
Technology  Laboratories,  Seattle,  U.S.A.  A  linear  array, 
sector  scanner  with  a  nominal  center  frequency  of  3.5  MHz 
and  a  field  of  vision  60°  was  used  as  the  transducer.  In 
the  cases  of  the  volunteer-I  and  the  patient,  the  focus  at 
transmission  was  set  at  a  depth  of  20mm  (approx,);  multi- 
focusing  had  been  used  for  the  volunteer-II,  with  focii  set 
at  depths  of  77,  100  and  125  mms  (approx.).  Dynamic  fo¬ 
cusing  was  applied  in  the  receiving  mode.  The  data  were 
sampled  at  12  MHz. 

IV.  Results 

A.  Cumulant  estimation  procedure 

Both  of  the  tissue  signatures  proposed  in  this  paper  rely 
on  computing  the  third-order  cumulants  from  digitized  RF 
scan  line  segments.  In  the  following  we  describe  the  pro¬ 
cedure  followed  to  estimate  them.  In  the  following,  we 
assume  y{k)  to  be  an  ergodic  process. 

(51)  Segment  the  i-th  data  record  yi{k)  into  K  records 
of  length  N  samples  each;  subtract  the  mean  of  each 
record  to  form  the  zero-mean  sub-segments  yij{k). 

(52)  Estimate  the  third-order  cumulants  Cy-.{r,p)  of 
each  sub-segment  yij{k),  j  =  0,2  •  •  -  Tf  —  1  according 
to: 

1 

Cy.i  i'Typ)  =  yij 

fc=0 

|r|<M,|/)|<M  (38) 

where  M  is  the  length  of  third-order  correlation  lags 
considered  in  the  computation. 
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(S3)  Average  the  cumulant  estimates  Cy.-{T^p)  over  all 
sub-segments  j  =  0, 1  ♦  •  •  A"  —  1  to  obtain  the  cumu¬ 
lant  estimate  Cy.  (r,  p)  corresponding  to  the  data  record 
yi{k)-. 

1 

CyX'r,p)  =  p)’ 

j=o 

B,  Estimating  the  periodicity 
B.l  Simulation  I 

In  this  experiment,  we  compare  the  proposed  method  to 
detect  periodicity  (see  eq.  (23)),  with  the  power  spectrum 
[10]  and  power  cepstrum  [26]  techniques.  Only  the  RPC 
component  of  the  RF  echo  was  considered.  We  simulated 
the  inter-scatterer  spacing  with  the  F-density  function,  for 

V  =  a/ p>  =  0.03,  0.15,  and,  0.35,  where  p  and  a  respectively 
denote  the  mean  and  the  standard  deviation  of  the  scatter- 
spacings.  The  nominal  mean  scatterer  spacing  was  p  = 
1.0mm. 

Cumulants  for  the  process  were  estimated  according  to 
the  procedure  (S1)-(S3)  with  the  following  parameters: 
M  =  100,  N  =  256,  A  =  4,  Fig.  1(a), (b)  and  (c)  show  the 
contour  maps  of  the  estimated  cumulants  corresponding  to 
=  0.03,0.15,0.35, 

Based  on  the  cumulant  estimates,  0  <  CJ  < 

TT,  was  computed  through  (23),  and  its  contour  map  is 
shown  in  Fig.  1  (d),(e)  and  (f)  for  v  =  0.03,0.15  and  0.35 
respectively. 

The  periodicity  in  data  is  evident  from  both  the  cumu¬ 
lants  and  in  p\{ri^uj).  However,  they  are  better  visualized 
in  Fig.l  (g),(h)  and  (i),  which  were  obtained  by  projecting 
the  absolute  values  of  Fig.  1  (d),  (e)  and  (f)  along  fre¬ 
quency  axes,  followed  by  smoothing.  The  dotted  lines  in 
Fig.  1  (g),(h)  and  (i)  correspond  to  the  true  mean  scatterer 
spacing,  1.00mm. 

Periodicity  can  also  be  estimated  from  the  w-axis  of 
p\{ri^u)\  Fig.  2  (a),  (b)  and  (c)  respectively  show 
0  ^  ^  <  tt,  corresponding  to  ri  =  0.0, 1.0  and 
2,0  mm  (see  Fig.l  (d),(e),(f)  and  eq.  (26)).  The  spacings 
between  adjacent  spectral  peaks,  which  is  given  by  27r/T 
according  to  (26),  can  be  obtained  from  Fig.  2  (a),  (b)  and 

(c)  and  periodicity  be  determined. 

Next  we  computed  the  power  spectrum  and  the  power 
cepstra  of  the  same  set  of  data.  Data  were  segmented  ex¬ 
actly  the  same  way  as  for  the  cumulant  estimates,  and  the 
average  autocorrelation  over  segments  obtained.  The  auto¬ 
correlation  estimates  were  smoothed  through  a  Blackman 
window,  before  computing  the  power  spectrum.  Fig.  3 
(a),  (b)  ,(c)  show  the  power  spectra  estimated  for  the  cases 

V  =  0.03,0.15  and  0.35.  Note  that  the  power  spectrum 
reveals  the  periodicity  only  for  the  case  v  =  0.03,  which 
corresponds  to  an  almost  periodic  scatter  spacings.  Based 
on  power  spectra,  the  power  cepstra  illustrated  in  Fig.  3 

(d) ,(e)  and  (f)  were  computed;  the  dotted  lines  in  these 
figure  corresponds  to  the  true  mean  scatter  spacing,  1.0 
mm.  The  power  cepstrum  shows  a  well  localized  peak  for 

=  0.03  and  a  reasonably  well  localized  peak  for  the  case 
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V  =  0.15.  However,  it  shows  a  large  spurious  peak  for  the 
case  V  =  0.35  which  overshadows  the  true  peak. 

We  ran  50  Monte  Carlo  runs  on  the  mean  scatterer- 
spacing,  based  on  the  power  cepstrum  and  the  proposed 
method,  corresponding  to  the  case  v  =  0.15,  and  the  mean- 
scatterer-spacing=  1.0mm.  Due  to  poor  performance,  the 
power  spectrum  method  was  excluded.  In  estimating  the 
periodicity  via  the  proposed  method,  the  first  peak  location 
was  obtained  as  the  location  corresponding  to  the  highest 
peak  in  the  region  0.4-2. 5  mm.  The  same  procedure  was 
used  in  the  case  of  the  power  cepstrum  method.  Table  I 
illustrates  the  results  we  obtained.  These  results  illustrate 
that  the  proposed  method  outperforms  the  power  cepstrum 
in  estimating  periodicity. 

Better  estimates  of  the  cepstrum  can  be  obtained  with 
parametric  methods  such  as  the  AR  approach  considered 
in  [26],  when  the  model  order  is  known  fairly  accurately. 
However,  the  estimation  of  the  model  order  is  not  always  a 
trivial  task.  In  this  paper,  our  use  of  non-parametric  tech¬ 
niques  to  compute  the  power  cepstrum  is  justified  because 
we  used  non-parametric  techniques  to  compute  higher- 
order  cumulants  too. 

B.2  Simulation  II 

In  this  simulation,  we  generated  the  RPC  component  as 
in  Simulation  I  and  added  a  DC  component  at  a  signal  to 
noise  ratio,  SNR  =  — 16  dB,  where,  following  [9]  the  SNR 
is  defined  as: 

where  Nr  is  the  number  of  resolvable  scatterers  in  the  data 
segment  considered.  If  the  normalization  term  Nr  is  not 
included,  the  SNR  measure  does  not  represent  the  true 
situation.  For  instance,  there  can  be  cases  with  only  a  few 
resolvable  scatterers  over  a  long  tissue  length,  where  the 
energy  per  scatterer  is  high  enough  so  that  the  periodicity 
can  be  determined  without  processing  the  signal,  but  still 
leading  to  a  misleadingly  low  SNR  measure. 

Fig.  4  illustrates  the  results  of  the  proposed  method  and 
the  power  cepstrum  one.  Fig.  4  (a),  (b),  (c)  show  |p3(Ti,a;)| 
projected  along  the  w-axis  for  cases  v  =  0.03,0.15,0.35. 
Fig.  4  (d),  (e),  (f)  show  the  power  cepstra  corresponding  to 

V  =  0.03,0.15,0.35.  Periodicities  obtained  from  Fig.  4  (a), 
(b),  (c)  are:  0.97,  1.00,  0.96  mm.  Periodicities  obtained 
from  Fig.  4  (d),  (f)  are:  0.96  and  0.88  mm.  The  result 
form  (e)  is  inconclusive;  there  are  several  spurious  peaks, 
two  major  ones  corresponding  to  periodicities  of  0.74  mm 
and  1.2  mm. 

Figs.  1,  4  and  Table  1  suggest  that  the  higher-order 
statistics  based  periodicity  estimation  technique  proposed 
in  this  paper  outperforms  both  the  power  spectrum  and 
the  power  cepstrum. 

A  problem  often  encountered  with  the  power  cepstrum  is 
that  it  tends  to  produce  spurious  peaks  in  the  presence  of 
noise.  Power  cepstrum  computations  are  heavily  affected 
by  the  DC  component  and  also  the  observation  noise  in 
RF  echoes.  Any  noise  escaped  through  the  averaging  and 
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windowing  operations  at  the  auto-correlation  estimation 
process,  may  contribute  a  large,  slowly  decaying  compo¬ 
nent  to  the  power  cepstrum,  thus  masking  small  peaks. 
Spurious  peaks  pose  a  serious  problem  to  the  periodicity 
detection/estimation  process.  Fig.  4  indicates  the  fact  that 
in  the  presense  of  the  DC,  the  problem  of  spurious  peaks 
gets  aggravated.  Since  HOS  suppresses  the  DC  compo¬ 
nent  and  measurement  noise,  the  proposed  technique  is  ro¬ 
bust  against  their  unwanted  effects. 

B.3  Clinical  data:  Normal  subjects: 

In  this  section,  we  present  the  periodicities  estimated 
from  clinical  images,  using  the  higher-order  techniques.  In 
all  clinical  images,  we  used  Af  =  80,  A  =  128  and  K  =  5. 
Each  segment  corresponded  to  a  length  of  8  mm  (approx.) 
in  tissue. 

Fig.  5  illustrates  results  corresponding  to  a  normal  liver 
(volunteer-I).  Contour  maps  of  0  <  w  <  tt  is 

shown  in  Fig.  5  (a),  and  its  smoothed  projection  along  the 
cj-axis  is  shown  in  Fig.  5  (b).  The  power  spectrum  and 
the  power  cepstrum  computed  from  the  same  data  using 
the  same  parameters  as  those  used  in  the  computation  of 
cumulants,  are  illustrated  in  Fig.  5  (c)  and  (d). 

The  period  estimated  from  the  cepstrum  and  the  pro¬ 
posed  method  are  respectively  1.28  mm,  and  1.3mm.  Pe¬ 
riodicity  is  much  more  evident  in  Fig.  5  (a),  (b),  which 
suggests  that  the  proposed  method  is  better  than  both  the 
power  cepstrum  and  the  power  spectrum.  From  the  power 
cepstrum/spectrum  it  is  impossible  to  see  repeated  peaks 
confirming  a  suspected  periodicity,  whereas  repeated  peaks 
are  clearly  seen  in  higher-order  estimates. 

Fig.  6(a)  shows  the  quadrature  demodulated,  loga¬ 
rithmically  compressed  liver  image  of  volunteer-I;  the  his¬ 
togram  of  estimated  periods  obtained  from  throughout  the 
image,  where  periodicity  could  be  detected,  is  shown  in 
Fig.  6(b).  The  average  value  of  the  period  evaluated  using 
the  proposed  method  over  80  estimations  was  1.03  mm,  at 
a  standard  deviation  of  0.26  mm.  Fig.  6(c)  and  (d)  show 
similar  results  for  the  volunteer  number  II,  where  the  av- 
erage=1.07  mm  and  the  standard  deviation  =  0.25  mm. 
In  these  images,  periodicities  were  estimated  in  a  depth  of 
approx.  35  mm  -  51  mm  from  the  transducer  surface. 

B.4  Clinical  data:  tumor  of  the  liver: 

In  this  section  we  report  the  results  of  our  attempts  to 
estimate  periodicities  from  focal  diseases  of  the  liver.  We 
obtained  a  clinical  image  of  a  patient  diagnosed  with  a 
hypo-echoic  tumor  of  the  liver  (see  section  3.2)  and  ap¬ 
plied  higher  order  periodicity  detection  schemes  proposed 
in  this  paper.  Fig.  7  (a)  shows  a  part  of  the  image  we  used. 
Periodicities  were  estimated  at  a  depth  of  48-65  mm  from 
the  transducer  surface.  The  focus  of  the  imaging  system  at 
transmission  had  been  set  at  about  20  mm  from  the  trans¬ 
ducer  surface;  dynamic  focusing  was  used  in  the  receive 
mode.  In  this  simulation,  we  used  M  —  80,  A  =  128  and 
K  =  5. 

Fig.  7(b)  shows  the  histograms  of  the  periodicities  es¬ 
timated  using  the  proposed  method  from  within  and  out- 
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side  the  tumor.  The  mean  and  the  standard  deviation  of 
the  scatterer  spacing  obtained  from  within  the  tumor  are 
respectively  1.07  and  0.42  mm.  For  regions  outside  the 
tumor,  corresponding  values  were  0.97  and  0.31  mm. 

Based  on  Fig.  7  we  conclude  that  there  is  no  signifi¬ 
cant  difference  in  periodicity  within  and  outside  the  tu¬ 
mor.  These  results  highlight  the  fact  that  there  are  im¬ 
portant  clinical  situations  where  the  periodicity  detection 
fails  to  delineate  normal  from  abnormal  tissue.  This  can 
happen  either  because  there  is  no  discernible  change  in  the 
periodicity  inside  and  outside  the  diseased  area,  or  there  is 
no  periodicity  associated  with  the  concerned  tissues  at  all, 
irrespective  of  normal  or  diseased. 

In  the  following,  we  investigate  the  performance  of  the 
color  of  the  tissue  response  as  a  tissue  signature  in  such 
situations. 

C.  Estimating  the  color  of  the  tissue  response 
C.l  Simulated  data 

We  simulated  two  realizations  of  digitized  RF  scan  line 
segments  corresponding  to  two  differently-correlated,  inter- 
scatterer  spacings.  In  the  first  case  an  RPC  component 
was  generated  with  t;  =  0.05.  The  simulated  RF  echo  was 
close  to  purely  periodic,  whose  cumulants  sequence  was 
only  slightly  affected  near  the  origin  by  the  slight  variations 
in  the  inter-scatterer  spacings.  This  served  as  our  first 
digitized  RF  scan  line  segment  under  study. 

In  the  second,  we  used  =  0.35  in  the  RPC  component 
resulting  in  a  digitized  RF  scan  line  segment  with  wider 
changes  in  the  inter-scatterer  spacings,  thus  contributing 
to  the  third  order  cumulant  near  the  origin,  i.e.  change 
the  statistical  color  of  the  RF  echo.  Another  component 
from  unresolvable  correlated  scatterers  were  simulated  with 
another  F  distribution,  withe  a/pi  =  0.8,  and  p  =  0.3 
mrUy  i.e.  inter-scatterer  spacing  distribution  was  almost 
non-periodic.  Correlations  were  further  introduced  to  the 
spacings  by  filtering  them  through  a  linear  time  invariant 
system  of  arbitrary  frequency  response.  The  echo  corre¬ 
sponding  to  this  scatterer  distribution  was  taken  as  the 
second  realization  of  the  digitized  RF  scan  line  segment. 

To  both  of  these  realizations,  a  DC  was  added  at  a 
SNR  =  defined  according  to  (40)  to  get  the  two 

digitized  RF  scan  line  segments  yi{k)  and  y2{k).  In  order 
to  simulate  the  effect  of  observation  noise,  we  added  white, 
Gaussian  noise  at  a  signal  to  noise  ratio  S/N  =  10d.B, 
where  the  ratio  S/N  was  defined  to  be: 

sllf }' ’  =  '■'  <«> 

where  yi{k)^  i  =  1,2  are  the  two  digitized  RF  scan  lines 
under  study  and  u;(^)  is  the  observation  noise. 

From  the  two  sets  of  digitized  RF  scan  line  segments 
yi{k)  and  y^ik)  so  generated,  we  estimated  the  third-order 
cumulant  sequences  c^^(Ti,r2)  and  according  to 

steps  51-54.  To  separate  the  non-periodic  parts  c^^(Ti,r2) 
and  C3^(ri,r2)  from  C3*(ri,r2),  i  =  1,2,  and  to  reduce 
the  estimation  variance  associated  with  the  cumulants,  we 
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applied  a  minimum  bispectrum  bias  supremum  window  of 
length  M  on  the  cumulant  sequence  obtained  through  steps 
51-54.  We  used  M  =  30,  A  =  128  and  K  =  8,  in  following 
51-54. 

Fig.  8(a)  shows  the  results  of  50  Monte  carlo  runs  for 
the  quantity  ui{k)  =  gi{k)  *  h{k)  estimated  through  the 
methods  of  section  2.2  ,  based  on  50  realizations  of  the  digi¬ 
tized  RF  scan  line  segments  corresponding  to  ?;  =  0.05,  i.e, 
yi{k).  The  average  value  over  50  runs  is  shown  in  a  solid 
line.  The  shaded  area  represents  the  average  ±  standard 
deviation.  The  dash-dotted  line  indicates  the  true  ultra¬ 
sonic  pulse  used  in  the  simulations.  Fig.  8(b)  shows  the 
corresponding  figure  obtained  with  50  realization  of  y2{k). 
It  can  be  seen  that  the  contributions  of  the  correlated  scat¬ 
terer  locations  affect  the  estimated  combined  kernel  U2{k). 
However,  from  Fig.  8  we  can  see  that  the  larger  quantity 
h{t)  masks  the  smaller  differences  between  Ui{k)  and  U2{k). 

We  formed  the  estimator  Sij,  i  —  l,2,j  1  according 

to  (32).  Fig.  8(c)  indicates  the  result  5ii  over  50  Monte 
Carlo  runs;  52 1  is  shown  in  Fig.  8(d).  The  average  estimate 
is  shown  in  a  solid  line,  while  the  shaded  area  represents 
average  ±  one  standard  deviation.  Fig.  8  suggests  that 
our  estimator  Sij,  i  —  l,2,j  =  1  is  capable  of  detecting 
the  color  difference  between  yi{k)  and  2/2 (Ar). 

C.2  Clinical  data 

To  investigate  the  performance  of  the  color  of  the  tissue 
response  as  a  signature  in  clinical  situations,  we  attempted 
estimating  the  color  of  the  tissue  response  measure  Sij  from 
images  of  the  liver  shown  in  Fig.  7(a).  In  estimating  the 
cumulants  we  used  iV  =  128,  M  =  25  and  K  ==  3,  where  the 
three  segments  corresponding  to  A"  =  3  came  form  three 
adjacent  A-lines  located  at  the  same  depth  of  the  image, 
i.e.  from  each  A-line  we  considered  data  segments  of  length 
128  samples  (8  mm  in  tissue  space). 

Fig.  9(a)  and  (b)  show  the  quantities  ui{t)  and  U2{t) 
estimated  respectively  from  within  and  outside  the  tumor. 
Solid  line  indicates  the  average  estimate,  while  the  shaded 
area  represents  average  ±  one  standard  deviation.  The 
estimates  5ii  and  52i  are  shown  in  Figs.  9(c)  and  (d). 

Fig.  9(c)  and  (d)  show  the  color  of  the  tissue  response 
estimated  from  within  and  outside  the  tumor.  In  form¬ 
ing  the  convolutional  inverse  Uj{k)~^  (see  32),  we  used 
the  average  of  20  combined  kernels  estimated  within  the 
tumor.  According  to  Fig.  9,  it  is  possible  to  delineate 
the  tumor  region  from  normal  tissue  using  the  estimate 
Sij,  i  =  1,2,  j  =  1.  The  closeness  of  the  estimator  to  a 
Dirac-delta  function  can  be  used  as  a  relative  tissue  signa¬ 
ture  in  the  present  case. 

Note  that  the  variations  in  the  estimations,  i.e.  the 
shaded  areas  in  the  figures,  are  not  entirely  due  to  statis¬ 
tical  estimation  errors.  It  has  been  shown  in  [1],  [2],  that 
using  the  same  higher-order  estimation  techniques,  under 
similar  conditions,  one  can  achieve  a  much  smaller  vari¬ 
ance.  A  major  reason  for  the  seemingly  larger  variance  is 
the  effects  of  natural  biological  variations  of  the  scatterer 
spacing  distributions  among  A-lines  across  the  image,  i.e. 
natural  differences  in  the  color  of  the  tissue  response  be- 
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tween  different  A-lines. 

V.  Discussion  and  Conclusions 

We  modeled  biological  tissue  as  a  collection  of  point 
scatterers  positioned  in  a  uniform  non-scattering  media. 
Based  on  the  higher  order  statistics  of  scatterer  spacings, 
we  derived  a  technique  to  estimate  mean-scatterer  spac¬ 
ing,  which  has  received  a  lot  of  attention  in  the  past  as 
a  potential  tissue  signature  for  organs  with  periodic  vari¬ 
ation  of  ultrasonic  properties.  We  also  proposed  a  higher- 
order-spectra  based  tissue  signature,  the  color  of  the  tissue 
response  to  supplement  situation  where  there  is  no  period¬ 
icity  or  periodicity  fails  to  delineate  normal  from  abnormal 
tissue. 

Using  simulated  data,  it  was  shown  that  the  periodicity 
can  be  better  visualized  and  detected  in  the  higher-order 
spectra  domain.  Since  zero-mean  Gaussian  processes  are 
suppressed  in  higher-order  domain,  the  proposed  method 
is  robust  against  additive  measurement  noise  unavoidable 
in  data  acquisition.  Microscopic,  diffuse  scatterers  of  the 
size  of  individual  cells  often  lead  to  i^F-echos  with  circu¬ 
lar  Gaussian  statistics.  Thus  the  proposed  method  is  in¬ 
sensitive,  unlike  power  spectra/cepstra  techniques,  to  the 
diffused  component,  because  of  the  property  that  HOS 
suppress  Gaussian  processes.  This  makes  the  higher-order 
spectra  domain  a  well  suited  place  to  estimate  periodicity 
associated  with  medical  ultrasound  data. 

The  periodicity  has  been  proposed  as  a  potential  tissue 
signature  in  detecting  liver  diffuse  diseases  such  as  cirrho¬ 
sis  and  chronic  hepatitis  [13],  [10].  Using  the  higher-order 
statistics  based  schemes  proposed  in  this  paper,  we  esti¬ 
mated  resolvable  periodicities  from  healthy  livers  and  a 
liver  with  a  hypo-echoic  tumor.  Although  we  were  able 
to  detect  and  estimate  periodicity  from  all  images,  we  ob¬ 
served  that  there  is  no  significant  difference  in  the  mean- 
scatter-spacing  estimated  from  within  and  outside  the  tu¬ 
mor  under  study. 

Detecting  and  estimating  periodicity  from  clinical  data  is 
not  always  a  simple  task.  In  methods  such  as  [10],  [13],  the 
presence  of  coherently  reflecting  blood  vessels,  for  instance, 
confuse  the  periodicity  estimation  process.  There  are  also 
a  lot  of  biological  variations  and  imaging-equipment  depen¬ 
dent  factors  such  as  the  focus  points  of  the  system,  beam 
width  etc.  which  compromise  the  definition  of  periodicity, 
and  its  estimation.  Tumors  are  highly  irregular,  complex 
structures  whose  organization  depends  on  the  type  of  the 
tumor  and  the  stage  of  growth.  There  can  be  regions  of 
necrosis,  and  regions  with  various  levels  of  perfusion,  served 
by  a  network  of  haphazardly  developed  blood  vessels  [23]. 
Even  when  the  organ  under  study  has  periodic  ultrasonic 
properties,  there  can  be  situations  where  there  is  no  char¬ 
acteristic  period  which  would  have  enabled  us  to  detect 
abnormalities.  It  is  plausible  to  expect  such  situations  in 
the  presence  of  a  malignant  tumor  growing  dendrites  into 
a  periodic  structure,  such  as  the  liver.  Methods  that  can 
only  estimate  periodicities  associated  with  resolvable  scat¬ 
terers  will  not  be  able  to  capture  the  true  picture  in  such 
cases. 
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The  higher-order  spectra  based  tissue  signature,  the 
color  of  the  tissue  response  proposed  in  this  paper,  com¬ 
pactly  summarizes  the  correlation  structure  existing  among 
the  scatterers  of  both  short-range  and  long-range  order.  It 
is  largely  independent  of  the  axial  pulse  echo  impulse  re¬ 
sponse  of  the  imaging  system  h{t)^  and  also,  does  not  re¬ 
quire  the  RF  echo  be  periodic.  Therefore,  the  color  of  the 
tissue  response  can  perform  as  a  tissue  signature  which  sup¬ 
plements  the  periodicity,  and  be  measured  even  when  the 
periodicity  is  not  defined.  Correlated,  unresolvable  scat¬ 
terers  can  be  expected  in  situations  such  as  the  tumor 
micro- vasculature.  Since  our  model  includes  the  correla¬ 
tions  among  unresolvable  scatterers,  the  effects  of  micro¬ 
vasculature  is  automatically  included. 

We  estimated  the  color  of  the  tissue  response  from  within 
and  outside  of  a  hypo-echoic  tumor  of  the  liver.  We  ob¬ 
served  different  colors  of  the  tissue  responses  for  the  two 
cases,  indicating  characteristic  scatter  correlation  struc¬ 
tures  inside  and  outside  the  tumor.  Currently,  investi¬ 
gations  are  under  way  to  devise  tumor  detection  schemes 
based  on  the  color  of  the  tissue  response,  as  observed  from 
RF  data  corresponding  to  traditional  B-scan  images. 
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TABLE  I 

Simulation- 1:  MONTE  CARLO  SIMULATIONS  FOR 
PERIODICITY  ESTIMATION  USING  POWER  CEPSTRUM  AND 
THE  PROPOSED  METHOD. 


Mean 

(mm) 

Variance 

(mm?) 

True  Period 

1.00 

0.000 

Power  cepstrum  method 

0.91 

0.026 

Proposed  method 

0.97 

0.005 

IEEE  TRANSACTIONS  ON  MEDICAL  IMAGING,  AUGUST  1996 


12 


Separation  (mm)  Separation  (mm)  Separation  (mm) 

Fig.  1.  Simulation  I:  (a),  (b),  (c)  Contour  maps  for  the  absolute  values  of  cumualnts  corresponding  to  v  =  0.03,0.15,0,35,  (d),  (e),  (f) 
contour  maps  of  |p3(tj  ,  u;)[,  0  <  o/  <  tt  for  cases  v  =  0.03,0.15,0.35,  (g),(h),(i)  smoothed,  projected  |p| (ri,  u;)[  along  the  u;-axis  for  cases 
V  =  0.03,0.15,0.35.  The  dotted  lines  in  (g),  (h)  and  (i)  indicate  the  true  periodicity,  1mm. 
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Fig.  2.  Simulation  I:  The  slices  of  the  |p3(ti,c<;)|  along  the  a>-axis  corresponding  to  (a)  ri  =  O.Omm,  (b)  r\  =  1.0  mm,  (c)  ti  =  2.0mm.  The 
spectral  peaks  are  separated  by  27r/T,  where  T  is  the  period. 


Frequency  (rad/S)  Frequency  (rad/S)  Frequency  (rad/S) 


Separation  (mm)  Separation  (mm)  Separation  (mm) 

Fig.  3.  Simulation-I:  (a),  (b),  (c)  The  power  spectrum  for  cases  v  =  0.03,0.15,0.35,  and,  (d),(e),(f)  the  power  cepstrum  for  cases 
V  =  0.03,0.15,0.35,  of  the  RPC-only  digitized  RF  scan  line  segments.  True  periodicity,  1mm,  is  indicated  in  (d),(e)  and  (f)  by  dotted 
lines. 
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024  024  024 

Separation  (mm)  Separation  (mm)  Separation  (mm) 

Fig.  4.  Simulation  II:  (a),  (b)  ,(c)  The  smoothed,  projected  along  the  a;-axis  for  v  =  0.03,0.15,0.35,  and,  (d),(e),(f)  the  power 

cepstrum  for  v  ~  0.03,0.15,0,35.  In  all  cases,  the  true  period  is  indicated  by  dotted  lines. 


Frequency  (rad/S)  Separation  (mm) 

Fig.  5.  Normal  liver:  (a)  The  contour  map  of  ct;)!,  (b)  the  smoothed,  projected  |p^(Ti,a;)|  along  the  cj-axis,  (c)  the  power  spectrum, 

and,  (d)  the  power  cepstrum  corresponding  to  the  liver  image  of  the  normal  subject,  Volunteer-!. 
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Fig.  6.  Normal  liver:  (a),(c)  B-mode  images  of  paxts  of  normal  liver  obtained  from  Volunteer-I  and  II,  (b),  (d)  the  normalized  histogram 
for  the  mean  scatterer  spacing  estimate  corresponding  to  (a)  and  (c).  The  number  of  estimates  involved  in  forming  the  histogreims  is  80. 
Average  and  the  standard  deviation  for  Volnnietr-I  are  1.03  and  0.26  mm.  Corresponding  numbers  for  Volunteer-II  are  1.07  and  0.25 
mm. 


Fig. 


7.  Tumor  of  the  Liver  (a)  The  B~mode  image  of  a  liver  with  a  hypo-echoic  tumor,  (b)  the  histograms  of  the  mean-scatterer  spacings 
estimated  from  inside  and  outside  the  tumor.  The  average  and  the  standard  deviation  estimated  inside  the  tumor  are  respectively  1.07 


and  0.42  mm,  while  the  corresponding  numbers  obtained  from  outside  the  tumor  are  0.97  and  0.31  mm. 
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Fig.  8.  Color  of  the  tissue  response:  (a)  The  combined  kernel  ui  (k)  corresponding  to  ~  0.05,  (b)  the  combined  kernel  U2  (k)  corresponding 
to  the  case  v  =  0.35  and  short-range,  correlated  scatterers.  The  true  pulse  used  in  the  simulation  is  shown  in  a  dash-dotted  line,  and 
the  average  estimate  evaluated  over  50  Monte  Carlo  runs  is  indicated  in  a  solid  line.  The  shaded  area  represents  average  ±  one  standard 
deviation,  (c)  The  estimator  5ii ,  and  (d)  S21 ;  solid  line  indicates  the  average  over  50  Monte  Carlo  runs,  and  the  shaded  area  represents 
the  average  ±  one  standard  deviation. 
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Fig.  9.  Tumor  of  the  liver:  Combined  kernels  estimated  form  (a)  inside,  and  (b)  outside  the  tumor;  (c)  the  color  of  the  tissue  response 
signature  Sn  estimated  inside  the  tumor  and  (d)  S21  estimated  outside  the  tumor.  In  all  cases  the  solid  line  indicates  the  average 
estimation,  and  the  shaded  area  represents  average  ±  one  standard  deviation. 


